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SECTION  1 
INTRODUCTION 


In  this  report,  we  describe  a 3-D  finite  difference  code,  DAVID, 
which  can  be  used  to  estimate  the  currents  and  voltages  induced  on  an 
arbitrarily  shaped  object,  located  over  a finite  conductivity  ground  plane 
if  desired,  when  illuminated  by  a plane  wave  gamma  source.  As  with  all 
3-D  codes,  the  spatial  resolution  that  can  be  obtained  is  severely  limited 
by  the  amount  of  computer  storage  available  and  the  amount  of  computer 
running  time  that  the  user  can  afford.  Also,  the  accuracy  of  the  physics 
must  be  compromised. 

DAVID  (and  DAVEJR)  was  developed  as  a research  code.  It  was  de- 
signed to  be  reasonably  accurate,  relatively  fast,  and  easy  to  understand 
and  modify.  It  is  intended  to  be  used  frequently  by  many  people.  Realiz- 
ing that  the  first  thing  a researcher  does  when  he  uses  a new  codr  is  "im- 
prove" it,  we  have  endowed  the  first  version  of  DAVID  with  only  the  most 
basic  physics  and  many  comment  cards.  The  time  and  spatial  steps  are  con- 
stant. Expanding  spatial  steps  are  impractical  with  arbitrarily  shaped 
objects  in  any  case.  We  try  to  make  up  for  the  loss  of  an  expanding  grid  by 
improving  the  outer  boundary  condition,  which  allows  it  to  be  closer  to  the 
object.  By  virtue  of  the  Cartesian  coordinate  system,  which  allows  us  to 
construct  an  object  in  a "building  block"  fashion,  the  field  equations  and 
electron  momentum  equations  are  in  their  simplest  form.  The  object  is  con- 
structed by  designating  certain  cells  of  the  grid  by  means  of  a flag.  This 
flag  causes  any  surface  tangential  electric  fields  and  normal  magnetic 
fields  to  be  set  equal  to  zero,  i.e.,  a perfect  conducting  boundary  condition. 
Any  fields  inside  the  object  are  also  zeroed. 
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Normally,  when  we  refer  to  DAVID  in  this  report,  we  will  also  be 
describing  DAVEJR.  The  only  differences  are  in  the  source  routine.  DAVID 
is  a particle  pushing  code,  so  that  self-consistent  effects  can  be  included, 
while  DAVEJR  uses  a prescribed  source  routine  with  a simple  modification  to 
approximate  electron  turning,  if  desired.  Without  the  particles  to  store, 
DAVEJR  can  be  made  to  run  much  faster  or  to  perform  higher  resolution  cal- 
culations. Each  code  is  useful  in  its  own  way.  The  physics  in  DAVID  is 
essentially  the  same  as  P0ST3D1  and  the  one-dimensional  phenomenology  code 
GLANC2’3. 

In  Section  2 we  will  discuss  the  basic  physics  that  is  involved  in 
the  close-in  coupling  calculation— including  those  aspects  which  we  do  not 
feel  can  be  handled  appropriately  in  DAVID.  In  Section  3 the  numerical  tech- 
niques are  displayed,  and  illustrative  calculations  are  shown  in  Section  4. 
Our  conclusions  and  recommendations  are  in  Section  5. 
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2.1  CLOSE-IN  PHENOMENA 
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The  phenomena  of  importance  to  close-in  EMP  coupling  are  those 
of  EMP  environment  prediction,  plus  boundary-layer  effects  and  surface 
electron  emission.  The  essential  physics  of  close-in  coupling  to  a vertical 
post  is  discussed  below. 

Consider  a vertical  cylindrical  post  protruding  from  a finitely 
conducting  ground,  as  shown  in  Figure  1.  Assume  that  the  line  of  sight 
from  the  post  to  a near-surface  nuclear  burst  makes  an  angle  0 with  the 
horizontal,  and  that  the  burst  is  sufficiently  removed  from  the  post  that 
the  gamma  wave  front  seems  planar. 
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As  the  gammas  from  the  burst  interact  with  the  air  they  produce 
primary  Compton  electrons  initially  moving  approximately  parallel  to  the 
gamma  flux.  At  points  well  above  the  ground  and  well  removed  from  the  post, 
pnly  electric  fields  are  initially  present.  As  time  progresses,  magnetic 
fields  are  generated  by  the  interaction  of  the  electric  fields  with  the 
boundaries,  and  the  primary  electrons  are  deflected  by  the  magnetic  fields. 
They  are  also  influenced  by  the  existing  electric  fields,  and  slowed  by  the 
effective  drag  force  due  to  ionizing  collisions  with  air  molecules. 

Ionizing  collisions  of  primary  electrons  with  air  molecules  create 
substantial  numbers  of  free  secondary  electrons  and  positive  ions.  Some  of 


Burst  (Far  Removed) 


Conducting 

Cylinder 


Problem  Boundary 

Figure  1.  Example  problem  geometry. 


the  secondary  electrons  attach  to  molecules  to  form  ions,  some  free 
electrons  recombine  with  positive  ions,  and  some  positive  and  negative  ions 
recombine.  Large  electric  fields  may  cause  electron  avalanching  in  the  air. 
The  rates  of  these  processes  are  all  distinct,  and  the  electron  attachment 
rate  to  0 depends  upon  the  electric  field  amplitude.  The  charged  species 
and  the  neutral  air  molecules  create  a collision-dominated  plasma.  Secondary 
electron  and  ionic  currents  may  thus  be  incorporated  into  Maxwell's 
equations  via  an  ohmic  conductivity.  Conductivity  is  calculated  using  the 
mobilities  of  electrons  and  ions,  where  electron  mobility  in  turn  depends  on 
the  amplitude  of  the  electric  field. 

Electromagnetic  fields  modify  the  primary  electron  trajectories 
as  well  as  the  electron  mobility  and  attachment  rates.  The  overall  problem 
is  clearly  nonlinear  and  must  be  solved  by  numerical  methods. 
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Self-consistent  treatment  of  the  primary  (Compton)  electron 
dynamics  is  the  key  problem  in  the  present  studies.  Many  of  the  phenomena 
outlined  above  can  be  included  in  three-dimensional  calculations  which  do 
not  treat  electron  dynamics  self-consistently.  However,  there  are  conditions 
where  even  the  initial  direction  of  current  flow  on  the  object  is  uncertain 
due  to  self-consistent  effects.  For  example,  consider  a gamma  flux  incident 
on  the  vertical  post  at  an  angle  of  about  30°  with  respect  to  the  horizontal. 
Magnetic  fields  due  to  conductors  tend  to  deflect  the  primary  electrons  away 
from  the  conductors.  The  ground  thus  tends  to  deflect  the  electrons  upward 
and  the  post  tends  to  deflect  them  downward.  Net  deflection  is  clearly 
uncertain,  as  is  the  initial  direction  of  current  on  the  post  which  would 
usually  oppose  the  vertical  current  in  the  air. 


Relative  responses  of  gamma-thick  and  gamma-thin  conductors  may 
also  be  greatly  modified  by  self-consistent  effects.  Self-consistent 
deflection  of  the  primary  electrons  may  reduce  the  charge  collection  by  the 
object.  For  gamma-thin  objects,  emitted  current  may  greatly  exceed  that 
collected,  while  for  gamma-thick  objects  the  total  charge  collection  may  be 
much  less  than  expected.  Emission  of  electrons  by  the  object  plays 
a pivotal  role  and  must  be  treated  as  accurately  as  possible.  We  do  not 
feel  that  a 3-D  code,  which  must  store  data  for  thousands  of  particles  being 
born  in  the  air  can  be  trusted  to  do  a reasonable  job  with  particle  emission 
from  an  object  as  well.  The  problem  is  compounded  by  the  lack  of  resolution 
in  the  spatial  grid.  Future  calculation  may  show  us  wrong,  but  as  far  as  the 
first  version  of  DAVID  is  concerned,  we  have  chosen  a different  approach  for  the 
case  of  a gamma-thin  object  in  air.  Instead  of  emission  specifically  from 
the  object's  surface,  we  allow  the  code  to  forget  that  the  object  is  there 
during  normal  particle  injection  and  movement  processes.  Unless  the  air  is 
very  thin,  this  is  a good  first  order  approximation  because  the  current  is 
almost  continuous  across  the  ooundary.  The  above  procedure  gives  us  a 
smoothly  varying  current  distribution  behind  the  object.  The  code  is  con- 
structed in  such  a way  that  it  would  be  almost  trivial  to  include  an  object 
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emission  scheme,  but  if  it  is  not  done  well,  the  unreal  current  distribu- 
tion in  the  layer  of  cells  behind  the  pole  can  give  rise  to  fields  which 
reflect  the  numerical  treatment  rather  than  the  physics. 

The  objects  used  in  DAVID  and  DAVEJR  can  be  either  gamma-thin 
or  gamma-thick,  i.e.,  they  can  be  either  completely  transparent  or  completely 
opaque.  If  we  were  limited  to  a single  simple  object,  e.g.,  a pole,  we 
could  easily  allow  for  an  object  with  a partial  shadow  (gamma  translucent). 
However,  when  one  part  of  an  object  can  shade  another  part  or  one  object 
can  shade  another,  the  logic  involved  with  partial  shadowing  can  take  up 
a significant  amount  of  computer  storage. 


DAVID  uses  a particle  treatment  of  the  Compton  electrons.  Particles 
representing  large  numbers  of  Compton  electrons,  are  injected  at  appropriate 
times  and  spatial  locations  within  the  problem  geometry;  weights  are  assigned 
to  the  particles  according  to  the  number  of  Compton  electrons  which  they 
represent.  All  of  the  particles  are  advanced  in  time,  using  the  Lorentz 
and  drag  forces  appropriate  to  each  individual  particle.  Based  upon  particle 
locations  and  velocities,  current  density  and  ionization  rate  are  calculated 
for  all  spatial  points  in  the  finite-difference  mesh.  The  air-ion  equations 
are  advanced  in  time  and  conductivity  is  calculated  for  each  point  in  the 
mesh.  This  is  done  using  the  existing  electric  amplitudes  at  that  point  to 
evaluate  the  field-dependent  mobility,  attachment  and  avalanche  parameters. 
Maxwell's  equations  are  then  advanced  in  time,  using  current  dersity  and 
conductivity  values  as  determined  above.  New  particles  are  injected  according 
to  the  time  and  spatial  distribution  of  the  gamma  flux.  The  process  is 
repeated  cyclically  until  the  desired  problem  time  is  reached. 


Because  of  computer  time  limitations,  the  present  state  of  the  art 
in  EMP  environment  calculations  cannot  be  realized  in  three-dimensional  close- 
in  coupling  calculations.  Treatments  of  gamma -ray  energy  spectra  and  initial 
angles  (and  angle-dependent  energy)  of  Compton  electron  ejection  lead  to 
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excessively  large  numbers  of  primary  particles  and  cause  unacceptably  long 
computation  time.  Instead,  monoenergetic  gammas  must  be  considered,  and 
the  primary  Compton  electrons  must  be  ejected  parallel  to  the  gamma  flux  (or 
normal  to  surfaces  in  the  case  of  surface  emission).  The  generation  of  a 
boundary  layer  near  the  surface  of  the  object,  wherein  positive  ions  ar.d 
electrons  separate  and  form  a gap  under  the  influence  of  the  normal 
electric  field,  has  not  been  treated  in  the  present  code.  Preliminary 
estimates  indicate  that:  (1)  the  high  capacitance  across  the  boundary  layer 

reduces  Its  electromagnetic  influence;  (2)  radiation  will  splash  electrons 
across  the  boundary,  reducing  its  influence  again;  and  (3)  the  contamination 
of  any  real  surface  will  affect  the  problem  in  such  a way  as  to  allow 
electron  charge  to  be  drawn  off  of  the  surface  more  easily  than  in  the  ideal 
case,  especially  with  the  added  influence  of  the  molecular  collisions  of  sea 
level  air. 


In  the  remainder  of  this  section,  we  present  and  discuss  the  actual 
equations  upon  which  DAVID  and  DAVEJR  are  based. 


2.2  FIELD  EQUATIONS 


The  field  equations  used  in  DAVID  are  Maxwell's  equations  in 
Cartesian  coordinates.  MKS  units  are  used  throughout,  with  the  exception 
that  the  magnetic  field  is  in  volts/meter,  i.e.,  the  quantity  ft  = ZQH  is 
carried,  where  Zq  is  the  impedance  of  free  space  (~ 120  ir  ohms).  For  a 
wave  propagating  in  free  space,  then,  the  electric  and  normalized  magnetic 
fields  would  be  equal.  The  use  of  ft  instead  of  H helps  in  studying  the 
physics  and  diagnosing  calculations.  In  conducting  regions,  we  have 
|e|  < |ft|.  In  order  to  keep  units  consistent,  the  current  density  must  also 
be  multiplied  by  ZQ. 

The  equations  used  in  DAVID  are,  in  vector  form 
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§ +ff  .c(Vx*-T)  , 


where 


?4, 


a » conductivity  (mho/m) 
e = permittivity  (farad/m) 

0 

c = speed  of  light  (3  x 10  m/sec) 
h = ZQH  (v/m) 

T = Z0J  (v/m2) 

H = magnetic  field  (amp/m) 

2 

J = driving  (Compton)  current  density  (amp/m  ) 
ZQ  = V p/e  (ohm) 
y = permeability  (henries/m) 


The  boundary  conditions  at  the  object  are  those  of  a perfect  con- 
ductor, i.e.,  the  tangential  E and  normal  H are  zero  (these  are  not  independent 
conditions) . We  assume  that  the  problem  has  mirror  symmetry  in  order  to 
decrease  the  number  of  grid  cells  required.  At  the  symmetry  plane,  the 
normal  E and  tangential  H fields  are  zero.  Two  types  of  outer  boundary 
condition  are  used:  the  perfect  conductor  (PC)  and  a fake  ambient  environ- 

ment (FAE)  condition.  The  FAE  condition  allows  one  to  move  the  outer  boundary 
much  closer  than  could  be  allowed  with  perfectly  conducting  walls  with  a given 
air  conductivity.  It  uses  some  of  the  field  characteristics  that  one  would 
see  if  the  object  were  not  present,  without  actually  calculating  those  fields. 
This  will  be  discussed  in  Section  2.  The  particular  treatment  used  in  DAVID 
can  be  improved  considerably,  but  has  been  shown  to  be  reasonably  successful 
even  in  its  primitive  form  (see  Section  4) . 


Even  though  a particular  boundary'  condition  may  not  cause  the 
currents  running  on  a single  pole  to  be  in  great  error,  it  does  change  the 
field  distribution  in  space  considerably.  The  pole  currents  do  not  change 
drastically,  at  early  times  at  least,  because  the  current  is  limited  in 
large  part  by  the  energy  stored  in  local  fields  located  very  close  to  the 
surface.  These  fields  determine  the  inductance  and  capacitance  of  the  pole. 
The  quasi-static  fields  are  not  greatly  affected  by  what  is  happening  far 
away.  However,  with  two  objects,  or  with  some  convoluted  object,  the  distri- 
bution of  the  fields  in  space  can  become  quite  important,  and  hence,  so  do 
the  boundary  conditions. 


2.3  MOMENTUM  EQUATIONS 


The  relativistic  Compton  electron  momentum  equation  in  our  system 
of  units  (h  = ZqH  = cB)  is 


■*  dE  £ 

dP  _ . v ft  e_  P 

dt  " " (E  c h)  - dRmf  p » 


where 


P = |?| 


/p2  + (me) 2 


dEe  e • 108o  (Ee  + 

dRmf  400  VEe  + 0'6) 


(Newtons)  , 


E = mc^  (Vl  + (P/mc)2  - 1)  (MeV)  , 
6 


*w9Ma «m -"■»  „.  ^ 


3 

p * air  density  (kg/m  ) 
m = electron  rest  mass 
me2  - 0.511 
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e = electron  charge  (1.6021  * 10  coulomb) 

The  drag  term  dEe/dRm^  is  obtained  by  differentiating  the  fitted 
mean  range-energy  relation*  and  converting  energy  to  joules: 

4.0E2 

Rmf  " p(EeVb.3)  (met<,rs>  ‘ (« 

Particles  are  advanced  once  each  time  step  using  previously  calcu- 
lated fields.  The  previous  value  of  P is  also  used  in  the  calculation  of 
the  drag  force  and  v/c. 

2.4  SELF-CONSISTENT  (PARTICLE)  SOURCES 

The  particle  motion  and  energy  loss  rate  must  be  converted  into 
Compton  currents  and  ionization  rates.  The  Compton  currents  go  directly 
into  the  field  calculation.  The  ionization  rate  is  the  driver  for  the  air 
chemistry  equation,  which  generates  the  electron  and  ion  densities  necessary 
for  the  conductivity  calculation.  The  conductivity,  in  turn,  is  used  in  the 
field  calculation. 

Each  particle  represents  the  number  of  electrons  formed  within  a 
cell  of  volume  dV  over  a period  of  time  dt.  The  electron  current  J is  then 
given  by 
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•>  2 
J * - eN  v (amp/m  ) , 


where  Nc  is  the  Compton  electron  density  (particles/nr)  and  v is  the  particle 
velocity  given  by  Equation  3a.  The  density,  Nc>  is  given  by 

Nc  = W/dV  (electrons/m3)  , (8) 

where  the  weight  W is  the  total  number  of  particles  and  is  given  by 

W = p p<t>dVdt  (electrons)  . (9) 

2 

Here,  p is  Compton  scattering  mass  attenuation  coefficient  (m  /kg),  p is 

c 2 
the  air  density,  and  <■)[>  is  the  average  photon  flux  (photons/m  - sec)  over 

the  time  interval  of  interest.  Actually,  since  DAVID  uses  the  same  cell  size 

everywhere,  the  dV  factor  is  not  necessary. 

In  order  to  save  storage,  particles  are  not  injected  at  every  cell 
so  an  averaging  scheme  is  necessary.  This  is  discussed  in  Section  3.4. 


The  mass  attenuation  coefficient  is  calculated  as  a function  of 
gamma  energy  from  the  scattering  cross  section  given  by  Evans4.  The  cross 
section  is 


a (cm2/elec)  = ur2Mj- £.n(l  + 2a)  + 


+'a)  (2a2-2a-l 
a2(l+2a)2 


3(l+2a)‘ 


where 


Tq  = classical  electron  radius 
= e2/mQc2  - 2.818  x 10"13  cm 
a = EY/ra0c2  * E /0.511 
E^  = gamma  energy  . 


The  advantage  to  using  the  mass  attenuation  coefficient  instead  of 
the  scattering  cross  section  is  that  it  is  essentially  independent  of  the 


•T 


material  or  its  physical  state4.  We  obtain  the  attenuation  coefficient  in 
the  following  way: 

2 2 

Pc(cm  /gm)  * e°cCcm  /electron)  x 7. 2 (electrons/air  atom)  * 

6.025  x io23  (atoms/mole)/14.4  (gm/mole)  . (11) 

The  conversion  to  MKS  units  is 

Pc(ra2/kg)  = 0.1yc(cm2/gm)  . (12) 


The  total  cross  section  (including  both  Compton  scattering  and 
absorption)  is 


e°  • -ol^r 


- 5 *»tlH ♦ h - ^2(lcm2/elei:) 


(13) 


The  absorption  cross  section  is  then 


a = o-a 

e a e e s 


(14) 


The  mass  absorption  coefficient,  p.  for  which  we  will  have  need  for  later, 

a 

is  calculated  in  the  same  way  as  p . 

v 

We  only  consider  Compton .processes  in  the  source  calculation. 

We  ignore  the  photoelectric  effect  and  pair  production.  This  will  be  reason- 
able if  we  confine  our  photon  energies  to  between  0.5  MeV  and  5 MeV. 


The  initial  direction  of  the  Compton  scattered  electron  is  taken 
to  be  parallel  to  the  direction  of  the  incident  gamma  rays.  The  electron 
is  given  an  energy  equal  to  the  average  energy  of  all  the  recoil  electrons 


T = E ( 0 / 0)  . 
av  y '•e  a e ' 


(15) 


The  electron  energy  is  about  j for  1.6  M«  ' gammas.  The  initial  electron 
momentum  is  then  given  by  the  inverse  of  Equation  5: 


AA 


^,v 


M -^  ,-*  ■.7  ”,''  ^ v * ; 


(16) 


The  ionization  rate  is  proportional  to  the  Compton  electron 
energy  loss  rate,  with  about  1 conduction  electron  being  created  for  each 
34  eV  lost  by  the  Compton  electron.  The  ionization  rate  necessary  for  the 
air  chemistry  equations  is 


Se(io„-P»irsM3-sec)  - Nc  ^ M°21  «.  lof  (17) 


mf 


3.4  x 10 


2.5  PRESCRIBED  (ANALYTICAL)  SOURCES 


Under  steady  state  conditions,  the  Compton  current  in  a medium 
which  is  homogeneous  over  the  electron  range  is  proportional  to  the  photon 
flux.  This  remains  true  with  a time  dependent  gamma  source  so  long  as  the 
electron  life  time  is  short,  so  that  equilibrium  is  maintained  and 
electromagnetic  fields  are  not  strong  enough  to  affect  electron  motion. 
These  two  conditions  generally  translate  into  high  material  density  and  low 
gamma  flux*. 


When  the  proper  conditions  are  present,  the  Compton  current  is 


equal  to 


(18) 


* In  the  case  of  a fast  exponentially  rising  pulse,  the  deviation  from  the 
steady  state  condition  manifests  itself  as  a simple  delay.  In  air,  this 
delay  is  several  nanoseconds.  In  the  ground,  it  is  negligible.  Therefore, 
properly  treated  prescribed  sources  would  have  the  ground  drivers  peaking 
before  the  air  Compton  currents.  This  delay  has  not  yet  been  built  into 
DAVEJR.  Its  absence  is  obviated  by  the  comparison  of  Figure  29. 
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where  t is  the  photon  flux  (photons/m  /sec) , i^  is  a unit  vector  oriented 
in  the  direction  of  the  photon  flux,  R ^ is  the  mean  electron  range,  given 
by  Equation  6,  and  is  the  gamma  mean  range.  DAVEJR  uses 

Ry  - (ycP)-1  (19) 

where  p is  the  air  density  and  uc  is  the  mass  Compton  collision  attenuation 
coefficient  as  described  in  Section  2.4. 


The  ratio  R^/E^R^  is  a fairly  constant  function  of  the  gamma 
energy,  E^5.  In  terms  of  this  ratio,  the  Compton  current  is 

J = • efY  rr , (20) 

Y EYRy 


where  the  gamma  energy  flux,  f is  given  by 

fY  ' V 1 


for  a monoenergetic  source.  Table  1 shows  the  ratio  for  various  values  of 
E^,.  It  is  a corrected  version  of  a table  used  in  Reference  5. 


A simple  two-piece  linear  fit  describes  the  ratio  well  over  the 
1-5  MeV  range*  The  fit  is  (see  Figure  2) 


Table  1.  Values  of  the  ratio  R f/E  R as  a func- 
tion of  gamma  energy. r ' y 


Ey(MeV) 


(1  /MeV ) 

Y Y 


0.0064 

0.0071 

0.0069 

0.0062 


This  curve  fit  c.nd  the  one  for  ionization  rate  are  not  used  in  DAVID/ 
DAVEJR,  but  are  shown  because  the  reader  may  prefer  to  use  this  alternate 
technique  in  a code  of  his  own. 


.J---.fi 


^^te-  * 


Rmf  ( 10  (5*°  + 1*4ey3  » 1 5 ey  - i-69  MeV 


E R 1 

Y Y (lO"3(7.9S  - 0.3SEy)  , 1.69  < E^  2 5 MeV 


The  ionization  rate  is  also  proportional  to  y.  In  DAVEJR,  we  use 
the  formula 


(Mi) . Map*  _ 


4 x 10 


where  u = mass  Compton  absorption  coefficient  (see  Section  2.4). 

a 

The  ratio  S /f  is  also  a slowly  varying  function  of  the  gamma 
e y 

energy  over  the  range  1 S E^£  5 MeV.  Table  2 shows  values  taken  from 
Reference  5 and  converted  to  MKS. 


A simple  curve-fit  describes  this  ratio.  It  is 


A/  ion.-2airs  \ = 115  E-o. 

f \m  - y - MeV/  fcy 


-0.275 


The  fit  and  data  points  are  shown  in  Figure  2. 

Prescribed  sources  are  used  for  underground  currents  in  both  DAVID 
and  DAVEJR.  The  parameters  are  the  same  as  in  the  air  (Equations  18  and  19). 


Table  2.  Values  of  the  ratio  Se/fY  as  a 
function  of  gamma  energy. 


Ey(MeV) 


S /f  ( i on-pai rs/m-y-MeV ) 

e Y 
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Figure  2.  Plot  of  ratios  useful  in  calculating 
prescribed  currents  and  ionization 
rates.  Curve-fits  are  layed  over  data 
points. 


The  photon  flux  is  attenuated  considerably  by  the  soil,  of  course  (see 
Section  2.6).  Provision  has  been  made  for  the  inclusion  of  a time  dependent 
soil  conductivity  caused  by  gamma -ray  enhancement.  There  is  a large  amount 


t 


■H 

J 

T 


of  uncertainty  as  to  what  the  dependence  of  the  enhanced  conductivity  on  the 
gamma  ray  flux  is.  At  this  time  the  codes  use  a function  which  makes  it 
proportional  to  y,  i.e.,  the  gamma  flux  (really  the  dose  rate,  but  the  two 
are  proportional  with  a single  photon  energy) . 

An  approximation  for  magnetic  field  turning  effects  on  the  Compton 


current  in  the  air  is  described  in  Appendix  A.  It  is  useful  only  in  the  | 
presence  of  small  fields,  when  the  electron  does  not  turn  too  far.  Electric  J 
fields  are  neglected.  The  latter  influence  could  be  included  with  no  dif-  . 1 
ficulty  using  the  same  principles.  At  this  time,  we  were  simply  looking  J 
for  the  first  order  effect  of  having  a new  component  of  current  introduced,  . H 
which  is  normal  to  the  radial  component.  Only  the  component  of  H which  | 
would  ordinarily  be  present  without  the  existence  of  an  object  (Hx)  is  4 
presently  considered.  The  other  components  could  easily  be  included  also.  1 
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2.6  THE  PHOTON  FLUX 

The  gamma  flux  at  a distance  r from  the  burst  is,  at  the  local 
time  t'  = t - r/c, 

,> 

•y(photons/m2-sec)  = f (t')  . (25) 

47TT  u 

The  function  F^Ct * ) is  the  gamma  time  history,  normalized  to  unit  area.  The 
time  histories  used  by  DAVID  and  DAVEJR  are  different  and  will  be  described 
later. 


It  is  assumed  that  the  variation  of  the  current  magnitude  with  r 
is  negligible  over  the  calculational  volume,  except  through  the  variation  of 
t'  and  except  for  attenuation  in  the  soil.  The  function  A(r),  in  the  air. 


A(r)  = exp(-MTPr)/EY  , (26) 

where 

= weapon  yield  (KT) 

Gy,  = gamma  efficiency 

K = 2.613  x io25  (MeV/KT) 

3 

p = air  density  (l.g/m  ) 

E^  = photon  energy  (MeV) 

2 

Up  * total  attenuation  coefficient  (m  /kg) 

r = range  over  which  gammas  are  attenuated 

The  function  A(r)  for  a point  in  the  ground  is  the  same,  except  that  it  is 
multiplied  by  an  extra  attenuation  term,  i.e., 

Vound  = A(r)  • exP('YgAr)  ’ (27) 


fSsmmir 
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where 


2 

p = ground  attenuation  coefficient  (m  /kg) 

5 

3 

p = ground  density  (kg/m  ) 

6 

At  = distance  traveled  through  ground 

The  codes  use  u - p , where  p is  the  same  as  that  of  air.  Note  that 
g c c 

multiple  scattering  terms  are  ignored.  These  are  important  in  determining 
the  current  distribution  at  small  incidence  angles,  i.e.,  nearly  horizontal. 


DAVEJR  uses  the  simplest  time  history  of  the  two  codes 

at 


F0(t)  - 


n O 

foe 


(a+b) (t-t _)  5 

1 + £e  0 

b 


(28) 


where  t^  is  the  time  of  peak  and  f^  normalizes  the  area  to  unity 


-at. 


f0  = ¥ (a+b)e 


(29) 


Since  FQ(t)  extends  back  to  t = it  is  necessary  to  adjust  tQ 
so  that  the  function  is  very  small  at  t = 0.  Given  the  ratio  R = F (0) / 

F (tQ) , an  approximation  for  tQ  is 


where 


(30) 


The  code  accepts  either  R or  t^  as  input. 
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DAVID  uses  the  more  complicated  four  piece  function  used  in  GLANC2. 
It  is  normalized  numerically  and  the  time  of  peak  is  an  input  parameter. 

The  function  is  shown  qualitatively  in  Figure  3.  The  parts  are  given  by 


iONH'.'?"'"""’"'*'  ’ 


. b (t-t_) 

Fj(t)  = Aje 


t < t 


1 » 


(31a) 


(lti) 


b„\  b_(t-t-) 
e z ^ 


r.2,  \ 3 1 

F0W  = A2  b2  (b2+bp(t-t2) 


, tj  < t < t 


2 * 


(31b) 


1 + zfr  e 
b3 


b2  Ktt-*?) 

1 + t ~e 

^3  ^ b3 

F0(t)  = A2  b2  (b^b3)  (t-t2) 
l + j—  e 
b3 


, t2  < t < t,  , 


(31c) 


Figure  3.  Qualitative  plot  of  the  y time  history 
used  in  program  DAVID. 
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In  order  to  make  this  function  the  same  as  the  simple  one  used  by 
DAVEJR,  set  the  following  equalities: 


'1 


b2  = b2 


b3  " b4 


In  addition,  set  t. 


tQ  and  maintain  < t2  < t3‘ 


Both  DAVID  and  DAVEJR  have  the  option  of  normalizing  the  y curve 
2 

to  a peak  flux  (MeV/m  -sec)  or  dose  rate  (rad/ sec) . If  the  yield  that  is 

4 

input  (YKT)  is  less  than  10  KT,  the  code  assumes  that  we  have  entered  a 

yield  in  KT,  a gamma  efficiency  (EPG) , and  an  attenuation  range  (ROB) . If 
4 19 

10  < YKT  < 10  , the  input  is  assumed  to  be  in  rads  (air)*  per  second. 

If  YKT  > lO1^,  the  input  is  assumed  to  be  in  units  of  MeV/m2-sec.  When  YKT 
is  in  rads/sec,  the  number  is  first  converted  to  MeV/m  -sec.  The  relation 
between  the  two  is 


(rad/ sec)  = 1.602  x 10"12  p fcm2/gm)F(MeV/m2-sec)  (31e) 

a 

where  p is  the  mass  Compton  absorption  coefficient  (see  Section  2.4).  The 
curve  is  normalized  to  unity  in  the  usual  manner,  but  the  coefficients  A^, 
A^,  and  A^  (DAVID)  or  f^  (DAVEJR)  are  then  renormalized  by  the  ratio  of  the 
desired  peak  "y  to  the  peak  y of  the  function  which  was  normalized  to  unity. 


2.7  AIR  CHEMISTRY 

In  order  to  calculate  air  conductivity,  the  densities  of  free 
(secondary)  electrons,  and  positive  and  negative  ions  must  be  known.  The 


* This  is  essentially  rad  (Si) /sec  for  our  purposes. 


- 


it 


28 


Sjjfca®?* 


V | 


treatment  of  these  quantities  is  identical  to  that  in  the  MRC  environment 
codes;  we  repeat  it  here  for  the  sake  of  completeness. 

Free  electrons  and  positive  ions  of  density  n and  n,  are  created 

3 e + 
at  a rate  Sg  (ion  pairs/m  -sec)  by  the  ionizing  collisions  of  piimary 

electrons  with  the  background  air.  Electrons  attach  to  0 2 with  rate  coef- 
ficient 8 forming  0”  ions,  of  density  n_.  Electrons  recombine  with  positive 
ions  with  rate  coefficient  a,  and  positive  and  negative  ions  recombine  with 
rate  coefficient  T.  Finally,  if  the  electric  field  strength  is  sufficiently 
high,  secondary  electrons  may  gain  sufficeint  energy  between  collisions  to 
ionize  air  molecules  in  subsequent  collisions,  creating  additional  secondary 
electrons  at  an  avalanche  rate  G. 
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We  assume  that  the  secondary  electrons  and  ions  everywhere  maintain 
local  charge  neutrality. 


n = n + n 
+ e 


(32) 


negative  ions  are 
dn 


ThT  = se  - (8  - G)ne  - anen+  , 
dn 

~ = 8ne  - Tn+n_  . 


The  assumption  of  local  charge  neutrality  is  clearly  not  satisfied 
in  a boundary  layer,  if  indeed  a boundary  layer  does  develop.  Further,  the 
usual  divergence  terms  are  not  present  in  Equations  33  and  34.  In  Equation 
33,  a term  V • nev^  would  normally  be  present.  However,  in  the  present 
problems,  secondary-electron  drift  velocities  are  such  that  the  distance 
of  characteristic  change  in  the  electron  drift  current  ngv£j  must  be  smaller 
than  10~4  meters  before  the  divergence  of  the  secondary  electron  current 
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becomes  comparable  to  attachment.  Clearly,  this  term  is  important  only  to 
the  formation  of  a boundary  layer  and  not  to  other  portions  of  the  problem 
where  quantities  change  characteristically  in  distances  on  the  order  of  0.1- 
1.  meters. 


Given  the  electron  and  ion  densities,  the  air  conductivity  is 
0 = e^ene  + Vn+  + nJ)  • (35) 

The  parameters  £S,  G,  a,  T,  yg  and  y^,  which  are  required  to  calcu- 
late the  air  conductivity,  are  obtained  from  curve  fits  to  existing  data 
performed  by  Longley,  Longmire,  Radasky  and  others.  The  parameters  are 
summarized  in  Table  3.  One  should  note  that  the  attachment  rate  $, 
avalanche  rate  G,  and  electron  mobility  y are  dependent  upon  local  electric 
field  amplitude,  relative  air  density  and  water  vapor  content  of  the  air. 

The  other  parameters  depend  upon  relative  air  density;  the  electron-ion 
recombination  rate  also  depends  upon  the  water  vapor  content  of  the  air. 
Field-dependent  parameters  for  dry  air  are  plotted  in  Figure  4. 


E/pr  (v/m) 

Figure  4a.  Electron  attachment  and  avalanche  coefficients 
for  dry  air. 
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SECTION  3 

THE  NUMERICAL  SOLUTION 


3.1  OVERVIEW  AND  PROBLEM  GEOMETRY 

Program  DAVID  is  written  as  a particle  moving  code  with  field  and 
air  chemistry  subroutines.  Since  the  particle  calculation  dominates  the 
running  time  requirements,  particle  information  is  stored  in  CDC  7600  small 
core  memory  (SCM),  while  the  field  and  ionization  arrays  are  stored  in  large 
core  memory  (LCM) . With  the  large  amount  of  storage  required  by  the  plotting 
packages’,  we  are  limited  to  about  4000  particles  in  SCM  at  one  time.  This 
usually  proves  more  than  adequate.  If  necessary,  we  could  buffer  groups  of 
4000  in  and  out  of  LCM.  The  CDC  6600  has  twice  as  much  small  core,  but  uses 
an  extended  core  memory  (ECM)  which  does  not  allow  random  access.  There 
would  be  the  usual  amount  of  difficulty  in  converting  from  one  machine  to 
the  other. 

Figure  5 shows  the  basic  program  flow.  Note  the  separation  of 
events  into  groups  for  which  numbers  appear  at  half-time  steps  and  integral 
time  steps.  Constant  time  steps  are  used  throughout  the  entire  calculational 
interval.-  The  program  flow  for  DAVEJR  is  the  same,  with  only  the  current 
calculation  being  different.  The  conductivity  calculation  is  shown  as  a 
separate  block.  In  the  codes,  it  is  included  as  part  of  the  electric  field 
calculation.  This  is  mostly  to  save  storage  (the  conductivity  is  not  stored 
and  must  be  recalculated  by  the  output  code) , but  it  is  also  convenient  to 
do  this  because  of  the  special  calculations  one  must  do  near  the  outer  boundary; 
the  fields  routine  was  already  designed  to  handle  the  boundaries  separately. 
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Figure  5.  General  flow  pattern  of  programs  DAVID 
and  DAVEOR. 


The  coordinate  system  used  in  the  programs  is  shown  in  Figure  6. 

Two  features  are  unusual.  The  first  is  that  the  vertical  coordinate  (y)  is 
pointed  downward  (which  must  be  remembered  when  looking  at  the  calculations 
of  vertical  field  and  current  components)  and  the  second  is  that  there  are 
two  separate  coordinate  systems:  one  in  the  air  and  one  in  the  ground.  In 

fact,  the  entire  solution  is  effectively  broken  into  two  regions  with  the 
techniques  being  different  in  each  case.  There  are  numerous  advantages  to 
having  the  coordinate  system  oriented  in  the  way  that  it  is,  and  these  may 
become  apparent  as  we  discuss  the  numerical  techniques. 

The  gamma  plane  wave  front  is  incident  from  the  left  and  upper 
sides  and  travels  in  the  +y  and  +z  directions.  Real  time  is  used  (as  opposed 
to  local  time)  and  t = 0 occurs  when  the  wave  front  reaches  the  y = 0,  z = 0 
line.  Note  that  all  events  occurring  along  on  x-coordinate  happen  at  the 
same  time. 
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Figure  6.  Coordinate  system  used  in  program  DAVID.  Note  the 

separate  coordinate  system  used  below  the  ground  plane. 


The  calculational  volume  is  divided  into  cells  with  the 
dimensions  Ax,  Ay,  and  Az.  Fields,  currents,  and  ionization  rates  are 
located  by  the  cell  in  which  they  are  assigned.  Figure  7 shows  how  the 
cells  are  arranged  within  the  grid.  Note  that  there  is  an  extra  layer  of 
cells  below  the  ground  surface.  These  extra  cells  contain  the  field 
components  necessary  for  specifying  the  field  boundary  conditions,  as  well 
as  acting  as  trash  cans  for  particles  leaving  the  grid.  In  the  future  they 
may  prove  useful  for  other  reasons,  such  as  improved  boundary  conditions  or 
for  holding  particles  which  represent  electrons  splashing  back  out  of  the  soil. 


Each  cell  has  six  field  components  associated  with  it,  as 
well  as  the  three  current  components  and  the  electron/ion  densities.  J 
and  the  densities  are  located  at  the  center  of  the  cell.  The  field  compon- 
ents are  located  along  the  sides,  as  shown  in  Figure  8.  The  electric  field 
components  are  centered  on  the  edges  to  which  they  are  parallel.  The  magnetic 
field  components  are  centered  on  the  sides  to  which  they  are  normal.  This 
has  obvious  advantages  when  specifying  boundary  conditions.  The  is  and  ft 

,35 


36 


I- 


v^\ 


Figure  8.  Unit  cell  used  in  program  DAVID.  Sides  are  numbered  to 
illustrate  the  grid  boundary  numbering  system. 


37 


M 


h 

5 


components  are  also  centered  with  respect  to  each  other  in  the  curl  sense. 

For  example,  the  calculation  of  E requires  the  z-component  of  V x H,  or 

3H  /3x  - 3H  /3y.  The  components  of  H that  are  needed  to  numerically  calcu- 
y x 

late  this  quantity  lie  on  opposite  sides  of  E so  that  the  curl  is  automatically 
centered.  The  currents  must  be  interpolated  between  cells  in  order  to  center 
them  on  the  electric  field  that  is  being  calculated.  The  same  is  true  of 
the  ion  densities  for  the  conductivity  calculation. 

Figure  8 also  shows  the  system  used  to  number  the  boundaries  of 
the  calculational  volume.  The  system  of  numbering  considers  three  types  of 
boundaries:  sides,  edges  (junction  of  two  sides),  and  corners  (junction  of 

three  sides).  Side  number  1 is  the  X = 0 plane;  side  number  2 is  the  X = XMAX 
plane,  etc.  Thus,  the  sixth  side  is  at  Z = ZMAX.  Edge  13  is  at  the  junction 
of  planes  1 and  3.  Corner  136  is  at  the  junction  of  edge  13  and  side  6,  or 
the  junction  of  planes  1,  3,  and  6.  This  may  all  seem  confusing  at  first, 
but  the  system  is  very  easy  to  remember  and  facilitates  working  with  the 
code— either  modifying  it  or  setting  up  a problem. 

Four  types  of  physical  boundary  conditions  are  used  in  DAVID.  They 
are  listed  in  Table  4.  The  variable  which  denotes  the  boundary  condition 
type  is  IBN  (NSIDE)  in  the  air  and  IBNG  (NSIDE)  in  the  ground.  NSIDE  denotes 
the  side  number,  as  described  above,  and  IBN  or  IBNG  take  the  values  1 
through  4,  depending  on  the  physical  condition  desired.  Not  all  sides  can 
assume  any  of  the  boundary  conditions.  Table  5 shows  the  types  of  boundary 
condition  that  each^side  was  allowed  to  assume  at  the  time  this  report  was 
written.  The  IBM  and  IBNG  arrays  are  read  as  input  data.  What  happens 
when  a forbidden  value  is  read  in  depends  upon  which  side  is  involved.  We 
will  now  briefly  describe  each  boundary  condition: 

1.  Perfect  Conductor.  Tangent  electric  fields  and  normal  mag- 
netic fields  are  zeroed.  Note  that  this  is  redundant,  since 
one  implies  the  other,  and  the  field  components  are  arranged 
within  the  cell  in  such  a way  that  zero  tangent  electric 
fields  automatically  yield  zero  normal  magnetic  fields, 
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Table  4.  Types  of  physical  boundary  conditions 
used  by  DAVID. 


IBN/IBNG 

Boundary  Condition 

1 

Perfect  Conductor 

2 

TM  Symmetry 

3 

Ground  Plane 

4 

Ambient  Air  or  Ground 

Table  5.  Types  of  physical  boundary  conditions  each 
side  is  presently  allowed  to  assume. 


Side 

Number 

Allowed  Values  Of 

IBN 

IBNG 

1 

2 

2 

2 

1,  4 

1,  4 

3 

1,  4 

3 

4 

1,  3,  4 

1,  4 

5 

1,  4 

1,  4 

6 

1,  4 

1,  4 
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2.  TM  Symmetry.  Also  known  as  "mirror  symmetry,"  this 
boundary  condition  effectively  doubles  the  size  of  a 
problem  wi+-h  a plane  of  symmetry.  In  the  future  it 
might  be  useful  to  give  plane  3 that  option  also.  In 
this  boundary  condition,  the  normal  derivatives  of  the 
tangent  electric  field  is  zero,  as  is  the  derivative  of 
the  transverse  magnetic  field.  In  practice,  it  is 
easier  to  set  the  TE  fields  equal  to  zero,  which  is  what 
is  done  in  DAVID.  In  order  to  do  *his,  the  boundary 
must  be  placed  through  the  center  of  the  cell  rather 
than  tangent  to  a cell  face.  Particles  must  also  be 
handled  in  a special  manner.  Whenever  a particle  tries 
to  cross  through  the  symmetry  plane,  its  normal  component 
of  velocity  is  reversed  in  order  to  represent  a particle 
coming  back  through  the  opposite  direction.  Special 
treatment  is  also  required  in  the  current  averaging 
scheme  near  the  symmetry  plane,  in  order  to  account  for 
the  contribution  of  electrons  on  the  opposite  side  of  the 
plane. 

3.  Ground  Plane.  Special  treatment  in  the  calculation  of 
the  tangent  electric  fields  at  the  air/ground  interface 
is  required  because  derivatives  of  the  magnetic  field 
across  the  boundary  are  required,  and  these  derivatives 
are  not  continuous.  We  decided  to  ignore  all  that  and 
just  take  the  derivatives  across  the  ground  plane. 

Special  handling  is  still  required  since  the  vertical 
grid  sizes  in  the  ground  are  different  from  the  air 
and  because  two  different  coordinate  systems  are  used. 

4.  Ambient  Air/Ground.  The  ultimate  objective  of  this 
boundary  condition  is  to  simulate  the  environment  that 
would  exist  in  the  absence  of  an  object.  The  best  way 
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to  do  this,  in  the  close-in  ground-burst  problem,  is  to 
use  a one-dimensional  calculation,  such  as  the  one  used 
in  GLANC2.  The  boundary  condition  is  important,  not  only 
for  specifying  the  electromagnetic  field,  but  for  the 
injection  of  particles  so  that  we  do  not  have  the  problem 
of  wasting  space  near  the  boundary  where  a sufficient  number 
of  particles  must  build  up  to  accurately  represent  the 
current.  At  this  time,  DAVID  does  not  handle  the  boundary 
condition  in  this  way.  It  uses  an  approximation  which 

exploits  several  characteristics  of  the  close-in  ambient 
field,  e.g.,  the  fact  that  the  TE  mode  fields  are  zero 
and  that  all  derivatives  in  the  x-direction  are  zero. 

There  is  the  additional  assumption  that  only  the  radial 
(direction  of  the  gamma  flux)  electric  field  exists  on  the 
top,  front,  and  back  boundaries  (sides  3,  5,  and  6).  That 
is  also  true  on  the  bottom  when  no  ground  is  present.  The 
x-component  of  the  magnetic  field  is  allowed  to  exist  on 
the  side  (side  2),  but  its  x-derivative  is  zero.  The 
boundary  condition  is  handled  somewhat  more  primitively 
in  the  ground. 


3.2  FIELD  EQUATIONS 

The  two  vector  field  equations  solved  by  DAVID  and  DAVEJR  in 
SUBROUTINE  FIELDS  are  (Section  2.2  Equations  1 and  2) 

H*  f E • c(V  * S - T) 

where  the  meaning  of  the  variables  is  given  in  Section  2.2.  For  the  pur- 
pose of  this  discussion,  we  will  use  capital  letters  instead  of  lower  case 


letters  for  the  normalized  magnetic  field  and  current,  i.e.,  h -+  H and 
j"  -*■  J.  The  reader  should  remember  that  the  quantities  are  normalized. 


Let  En^i»  k)  be  the  uth  component  (x,  y,  z)  of  the  electric 
field  evaluated  at  the  n time  step  and  in  the  cell  whose  x,  y,  and  z 
coordinates  (cell  center)  are  [(i-l/2)'Ax,  (j-l/2)Ay,  (k-l/2)Az].  The 
notation  for  the  H-field  is  similar.  Remember  that  £ and  ft  are  evaluated 
half  a time  step  apart.  Thus  e"  is  half  a time  step  (At/2)  later  than 
Hj.  Similarly,  the  positions  of  the  field  components  are  different,  even 
though  they  are  designated  by  the  same  i,  j,  k indices.  In  the  code,  the 
solution  for  the  electric  field  occurs  earlier  in  the  loop  than  the 
magnetic  field  z because  H is  assumed  to  be  zero  during  the  first  pass. 


The  magnetic  field  equations  are  center  differenced.  They  are 


H"+1(i,j,k)  = H^(i,j,k)  + cAt 


~Ey(i.J> 

- 


k+l)  - E“(i,  j,k) 


E"(i,j+l,k)  - E“(i,j 


n+l  n fci+l,}, 

H“  1(i,j ,k)  = Hj(i,j,k)  + cAt  — 


i.j.k)  ' 

y 


k)  - E^(i,j,k) 


E^(i,j ,k+l)  - E"(i,j 


9 


n+l  n Wd.M.k) 

H”  X(i,j,k)  = Hj(i,j,k)  + cAt  — 


E"(i+l,j,k)  - E"(i,j 


i»j,k) 


(43) 


£“Ci,3,k)  = 


Hy+1(i,j,k)  - Hy+1(i-l, j ,k) 
Ax 


H"+1(i,j,k)  - H^fi.j-l.k) 


Ay 


- J 


n+1 


The  current  components  are  calculated  at  the  cell  centers  and  must 
be  interpolated  to  find  their  values  at  the  position  of  each  electric  field 
component.  Similarly,  the  conductivity  must  be  calculated  at  each  E-field 
position,  so  that  the  electron  and  negative  ion  densities  must  be  inter- 
polated. We  illustrate  the  procedure  for  the  case  of  the  current  components, 
but  the  scheme  is  the  same  for  the  ions. 


4nt  = \ [Jx(i,j,k)  + Jx(i,j~l,k-1)]  , 
Jynt  tyM.k)  + Jy(i-l,j,k-l)]  , 


= \ [Jz(i,3,k)  + Jz(i-l,j-l,k)]  .' 


(44) 

(45) 

(46) 


The  conductivity  calculation,  which  takes  place  in  the  FIELDS 
routine,  is  discussed  in  Section  3.7. 


Note  the  form  in  which  Equation  40  was  written.  The  term  in  brackets 


is 


-anAt 

u 


n«4. 

a At 
u 


should  mathematically  go  to  the  limit  of  unity  as  a^At  goes  to  zero. 

This  occurs  for  very  low  conductivity.  In  order  to  allow  the  code  to  operate 

in  the  limit  of  zero  conductivity  (free  space  propagation)  we  add  a small 

—4  n 

constant  (10  ) to  the  dimensionless  quantity  auAt.  This  gives  us  a 


44 


-f  _ "jisr  1 <T— 


K-v»^«.'-.f--»  v-.-~, , . _. 


s 

& 


f 


&*' 


if 


fer 

&■ 

®; 


ss. 

| 

t" 

-’S’*' 

K 


I 


numerical  result  sufficiently  close  to  unity  and  does  not  allow  the  computer 
to  try  to  divide  by  zero. 

The  field  equations  are  differenced  over  the  entire  mesh,  regard- 
less of  the  presence  of  the  object.  After  each  field  is  differenced,  we 
loop  back  through  the  mesh  and  erase  the  fields  within  the  cells  correspond- 
ing to  the  body  and  the  tangential  electric  or  normal  magnetic  field  components 
on  the  surface. 

Special  interpolations  are  required  to  obtain  the  currents  and 
ion  densities  for  the  calculation  of  fields  near  a boundary  because  one 
cannot  interpolate  through  it.  In  order  to  avoid  having  many  "IF"  checks, 
DAVID'S  field  subroutine  calculates  the  fields  along  each  side,  edge,  and 
comer  explicitly.  This  requires  a lot  more  programming,  but  by  making  the 
decisions  beforehand,  instead  of  letting  the  computer  do  it,  a significant 
amount  of  computer  time  can  be  saved. 

As  will  be  seen  in  the  next  section,  particles  and  current  averaging 
near  boundaries  must  also  be  treated  in  special  ways.  For  all  these  reasons,  we 
have  assigned  a flag  to  each  cell  in  the  mesh.  This  array  is  called  IBOD 
for  the  air  cells  and  IBODG  for  the  ground  cells.  Each  element  of  the  array 
has  the  value  0 through  27.  A cell  whose  IBOD  (or  IBODG)  is  0 is  not  in  an 
object  or  next  to  a boundary.  A cell  whose  array  element  is  numbered  1 through 
6 has  one  side  near  a boundary  (the  number  tells  which  boundary  it  is).  Cells 
numbered  7 through  18  have  two  sides  on  a boundary  (located  on  an  edge) . 

Cells  numbered  19  through  26  have  three  sides  adjacent  to  boundaries  (located 
in  corners).  Finally,  cells  with  IBOD  or  IBODG  elements  equal  to  27  are 
part  of  the  object.  Thus,  an  object  is  constructed  by  setting  IBOD  and  IBODG 
equal  to  27  for  all  the  cells  which  one  wishes  to  use  to  describe  it.  The 
object  takes  on  the  appearance  of  a model  built  out  of  blocks.  The  IBOD 
numbering  system  is  summarized  in  Table  6. 
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Table  6.  IBOD/IBODG  numbering  system. 


IBOD/IBODG 

Cell  Location 

0 

Away  from  boundary 
and  object 

1 - 6 

Next  to  side 

7 - 18 

In  two-sided  edge 

19  - 26 

In  three-sided  corner 

27 

In  object 

A flow  chart  showing  the  electric  field  calculation  is  shown  in 
Figure  9.  Note  that  the  conductivity  is  calculated  during  the  electric 
field  calculation  (SUBROUTINE  FIELDS)  and  is  not  stored.  The  time  output 
code,  DAVEOUT,  contains  a routine  to  recalculate  the  conductivity  from  the 
ion  densities  and  the  electric  field.  Thus,  the  outputted  conductivity  is 
somewhat  different  than  the  value  actually  used  in  the  E-field  calculation. 

With  the  exception  of  a relatively  thin  region  near  the  surface, 
the  conductivity  of  the  ground  remains  constant.  It  is  only  within  that 
same  layer  that  currents  of  significant  magnitude  exist.  In  order  to 
speed  the  calculation  somewhat,  the  ground  is  divided  into  two  regions.  The 
one  near  the  surface  uses  a conductivity  array  and  a current  array  which 
change  in  time.  This  layer  is  NRY  cells  deep.  The  array  is  two-dimensional; 
no  x-direction  variation  is  allowed.  Prescribed  sources  are  used  and  the 
gamma  flux  is  attenuated  exponentially  with  slant  range.  The  field  equations 
for  the  remainder  of  the  grid  do  not  include  a time  dependent  radiation 
enhanced  conductivity  or  a Compton  current  term.  The  subroutine  which 
calculates  the  ground  electric  fields  is  called  GELEC  and  the  magnetic  field 
routine  is  called  GMAG.  The  names  of  most  variables  in  the  ground,  including 
the  coordinates,  are  the  same  as  in  the  air,  except  that  the  letter  "GM  is 
added . 


Figure  9.  Flow  chart  of  the  calculation  of  electric 

fields  in  the  air  as  performed  by  subroutine 
fields. 
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In  the  future,  the  code  should  be  modified  to  allow  smaller  cells 
in  the  radiation  deposition  layers,  than  in  the  remainder  of  the  ground. 

A description  of  the  ambient  boundary  condition,  as  presently 
implemented,  is  in  order.  See  Section  3.1  for  definitions  and  numbering 
systems.  We  start  with  the  air  calculation. 

Plane  1,  the  TM  symmetry  plane (X  = 0),  runs  through  the  center  of 
the  first  cells  (I  = 1).  Thus,  the  field  components  EX(1,J,K),  HY(1,J,K), 
and  HZ(1,J,K)  are  located  on  the  plane  and  are  equal  to  zero.  The  TM 
fields  in  these  cells  are  set  equal  to  the  fields  in  the  next  cell  in  the 
x-direction,  i.e.,  EY(1,J,K)  = EY(2,J,K),  EZ(1,J,K)  = EZ(2,J,K),  and  HX(1,J,K) 
HX(2,J,K) . 

The  fields  on  plane  2,  X = XMAX  (I  = NX),  are  calculated  using  the 
full  field  equations  for  EY(NX,J,K)  and  EZ(NX,J,K),  except  that  partial 
derivatives  with  respect  to  X are  assumed  to  be  zero. 


On  the  top  plane,  plane  3 (Y  = 0 and  J = 1) , EX(I,!,K)  = 0,  HY(I,1,K) 
0 and  EZ(I,1,K)  are  calculated  as  components  of  the  radial  electric  field 
(curl  H terms  zero) . We  impose  the  condition  that  EZ  has  no  variation  in 
the  x-direction  by  calculating  EZ(2,1,K)  and  using  this  value  for  I 5 2. 


The  calculation  at  plane  4 (Y  = YAMAX,  J = NY)  is  the  same,  if 
the  ambient  boundary  condition  is  used.  An  air/ground  interface  may  also 
be  used,  in  which  case  we  difference  across  the  ground  plane  using  magnetic 
fields  from  the  ground  field  array. 

On  planes  5 and  6 (Z  = 0,  ZMAX  and  K = 1,  NZ)  EX  is  identically 
zero  and  EY(I,J,1)  and  EY(I,J,NZ)  are  calculated  components  of  the  radial 
electric  field.  As  on  plane  3,  the  x-variation  is  set  equal  to  zero  by 
calculating  EY(2,J,1  or  NZ)  and  using  this  value  for  I 5 2. 
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The  ambient  boundary  conditions  in  the  ground  are  somewhat  dif- 
ferent because  there  is  no  dominant  radial  electric  field.  The  symmetry 
plane  (plane  1)  is  treated  exactly  as  it  is  in  the  air.  The  EZG(NX,J,K) 
and  EYG(NX,J,K)  fields  on  plane  2 are  set  equal  to  those  just  inside  the 
boundary  (I  = NX-1).  The  fields  on  plane  3,  the  air/ground  interface  are 
set  equal  to  those  in  uic  air,  i.e.,  EXG(I,1,K)  = EX(I,NY,K),  EZG(I,1,K)  = 
EZ(I,NY,K),  HYG(I,1,K)  = HY (I ,NY,K) . 

On  planes  4,  5,  and  6,  a constant  of  proportionality  is  used  to 
determine  the  fields  on  the  boundary  from  values  inside  the  mesh.  The  con- 
stant, called  GUESS,  is  currently  set  equal  to  0.9.  If  it  were  set  to  0., 
it  would  be  equivalent  to  an  infinite  conductor  boundary  condition.  The  use 
of  GUESS  suppliments  the  requirements  that  partials  with  respect  to  x are 
zero  and  that  EXG,  HYG,  and  HZG  are  zero  on  the  boundary.  For  example, 
on  plane  5 we  have  EXG(I,J,1)  = 0.  For  I = 2,  we  set  EYG(2,J,1)  «*  GUESS* 
EYG(2,J,2)  and  EZG(2,J,1)  = GUESS*EZG(2 ,J ,2) . Then,  to  maintain  zero 
derivative  in  the  x-direction,  we  set  EYG(I,J,1)  = EYG(2,J,1)  and  EZG(I,J,1) 
EZG(2,J,1)  for  I > 2. 

3.3  MOMENTUM  EQUATIONS  AND  AVERAGING  TECHNIQUE 

The  particle  calculation  is  the  most  complicated  of  all  those 
performed  by  DAVID  and  forms  the  main  part  of  the  code.  The  flow  of  the 
computation  is  shown  in  Figure  10. 

The  process  starts  with  the  injection  of  new  particles  in  alternat- 
ing cells,  like  a 3-D  checkerboard.  The  user  can  choose  whether  he  wants 
to  inject  each  time  step  or  to  skip  one  or  more.  The  input  variable  INCINJ 
does  this  (INCINJ  = 1 injects  every  time  step).  In  order  to  help  smooth 
the  current  distribution,  the  cells  in  which  injection  ccci.rs  alternate  each 
injection  time.  For  example,  if  injection  occurred  in  the  1st,  3rc*,  5^, 
etc.,  cells  in  the  x-direction  during  the  first  time,  the  2nc*,  4t^,  6t*1, 
etc.,  cells  would  be  injected  during  the  second  time. 
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TO  INJECT? 


INJECT  PARTICLES  IN  ALTER- 
NATE CELLS  (ALSO  ALTERNATING 
IN  TIME).  STORE  PARAMETERS 
IN  ENDLESS  LOOP.  DU  NOT  EMIT 
IN  OBJECT  OR  SHADOW  IF  OBJECT 
IS  THICK 


MAXW 

EOUA" 

_VE 

ELL'S 

HONS 

ENTER 

LOOP 


I 


RECALL  OLD  MOMENTA, 
POSITION,  AND  PAR- 
TICI  E WEIGHT 


T 


S5 


— WJl 

AIR  CHEMISTRY 
EQUATIONS 


LOOP  THROUGH 
PARTICLE  STACK 


8 


) 


EXIT 

LOOP 


RETURN 


COMPUTE  PARTICLE 
ENERGY  AND  DRAG 
FQPCE 


Figure  10.  Flow  chart  showing  the  particle 
motion  and  current/ionization 
rate  calculation. 
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Figure  10  (continued). 
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Figure  10  (continued) 


PY (y-momentum) , PZ (z -momentum) , X(x-position) , Y (y-position) , Z(z-position) , 
and  W(particle  weight).  The  initial  position  of  the  particle  is  the  center 
of  the  cell  in  which  it  is  injected.  The  initial  kinetic  energy  of  the 
particle  is  given  by  Equation  15  and  its  momentum  is  then  calculated  by 
Equation  16.  The  particle  weight  is  given  by  Equation  9. 

The  PRW  array  is  circular.  When  the  last  particle  is  entered 
(4000  are  currently  allowed),  the  code  writes  over  the  first  ones  until  all 
new  particles  are  entered.  The  first  ones  are  usually  the  older  ones. 


When  particles  die,  the  surviving  ones  are  pushed  up  the  stack  so  that  new 
ones  can  be  entered  below.  It  is  only  when  new  ones  can  no  longer  be  entered 
there  that  we  replace  the  old  ones  at  the  top.  It  is  not  foolproof,  but  it 
seems  to  be  the  best  of  the  simple  schemes. 

When  the  object  has  been  designated  as  gamma  thin  (IBLACK  = 0), 
particles  are  injected  everywhere,  including  the  cells  inside  the  pole. 

They  are  also  allowed  to  travel  through  the  pole  when  born  outside  of  it. 

If  the  object  is  gamma  thick  (IBLACK  = 1)  no  particles  can  be  born 
in  it  or  in  its  shadow.  An  array  called  SHADO  (I,J,K)  in  the  air  and  SHADG 


(I,J,K)  in  the  ground  contain  a shadow  factor  for  each  cell.  The  shadow  % 
factors  are  calculated  in  SUBROUTINE  SETUP  from  the  object  description  and  J 
the  incidence  angle  of  the  gamma  rays.  If  the  center  of  the  cell  is  within  % 


the  shadow,  the  entire  cell  is  considered  to  be  within  the  shadow. 

Ordinarily,  the  shadow  factor  either  has  the  value  0 or  1.  The  value  0 means 
t!.e  cell  is  completely  shadowed  and  the  value  1 means  that  it  is  not.  When 
the  ambient  boundary  condition  is  used,  we  gradually  increase  the  shadow  factor 
for  cells  near  the  boundary  which  would  otherwise  be  shadowed.  If  we  did 
not,  there  would  be  a contradiction,  since  a shadow  could  not  exist  in  the 
ambient  environment.  The  integer  KSTP  determines  the  number  of  cells  over 


next  to  the  wall  has  the  value  0.666* ••  and  the  next  one  has  the  value  0.333 
and  the  third  one  is  zero.  When  the  shadow  factors  value  is  between  0 and  1, 
it  is  used  to  multiply  the  particle  weight. 

Particles  are  not  injected  until  the  incident  gamma  wave  front 
passes  through  the  cell  center.  Since  the  t waveform  rises  exponentially 
from  t = - 00 , we  include  all  of  the  charge  that  would  be  generated  over  that 
time  period  in  the  weight  of  the  first  particle. 

After  injection,  the  particle  momentum  calculation  begins.  Currents 
are  not  calculated  until  the  particle  has  been  affected  by  drag  and  fields. 
Starting  with  the  first  particle  in  the  stack,  the  old  electron  energy  (EE)  is 
calculated  using  Equation  5.  This  allows  the  calculation  of  the  drag  force 
(Equations  4 and  3).  The  indices  of  the  cell  (ICELL,  JCELL,  and  KCELL)  are 
computed  from  the  x,,y,  and  z coordinates.  The  electric  and  magnetic  fields 
at  the  center  of  the  cell  are  computed  by  interpolation.  Then,  the  momentum 
Equation  3 is  integrated,  using  the  old  value  of  momentum  on  the  right-hand 
side  in  the  drag  term. 


The  equations  for  the  three  components  of  momentum  can  be  put  in 
the  form  of  Equation  39  and  the  exponential  form  of  solution  used.  This  is 
useful  at  low  altitudes  where  the  drag  term  is  significant.  As  an  example, 
we  write  the  solution  for  the  x-component  of  momentum 

-A  'A 

PXn+1  = PXne  C - eFxAt — ) > (47) 

where  we  are  calculating  PX  at  the  (n+1)1"*1  time  step  and  where 
e = electron  charge 
At  = time  increment 


. dE  At 
c " dR  P 
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P = old  momentum  magnitude 


F = EX  + 
x 


PY-HZ  - PZ-HY 
\fp2  + (mc)2 


HZ,HY  = interpolated  magnetic  field  (z  and  y components) 
(volts/m) 

EX  = interpolated  electric  field  (x-component) 
m = electron  rest  mass 
c = speed  of  light. 

When  the  new  momentum  components  are  computed,  the  new  velocities 
can  be  found  using  Equation  3a,  and  then  the  new  position  can  also  be  computed. 

If  the  particle  has  penetrated  the  symmetry  plane,  its  x-momentum 
is  reversed  and  it  is  placed  an  equal  distance  on  this  side  of  the  plane 
(plane  1) . 

The  charge  density  (QDEN)  is  calculated  from  W,  and  the  three 
components  of  current  are  determined  as  the  products  of  it  and  the  velocity 
components.  The  new  energy  is  then  computed  and  the  ionization  rate  (S) 
figured  from  Equation  17. 


DAVID  then  determines  whether  the  old  position  of  the  particle 
was  near  a boundary  wall.  This  is  done  by  inspecting  IBOD  (ICELL,  JCELL, 
KCELL).  If  the  particle  was  in  the  cell  near  a boundary  before  moving  this 
time  (1  S IBOD  5 26)  it  may  have  penetrated  and  must  be  killed  after  its 
current  has  been  spread  in  a manner  appropriate  to  that  particular  boundary. 
Most  particles  will  not  be  in  one  of  those  cells  and  one  should  not 
waste  time  by  dropping  out  of  the  loop  or  doing  special  checks.  That  is 
why  the  IBOD  array  is  used.  If  IBOD  = 0 or  27  the  particle  was  clear  or 
in  a pole  away  from  the  boundary  (transparent  pole  only,  since  the  particle 
would  already  have  been  killed  if  it  had  run  into  a thick  pole),  if  IBOD 
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does  not  have  these  values,  we  drop  out  of  the  loop.  A 26  position  computed 
GO  TO  uses  the  value  of  IBOD  (ICELL,  JCELL,  KCELL)  to  send  the  pointer  to 
the  part  of  the  program  that  will  handle  the  particle  at  the  particular  side, 
edge,  or  corner.  Ordinarily,  the  code  will  qnter  directly  into  the  current 
smearing  phase. 

Because  constant  spatial  step  sizes  are  used,  a rather  unsophisti- 
cated averaging  scheme  can  be  used  to  smear  the  current  and  ionization  rate 
in  space;  one  does  not  have  to  remember  what  the  volume  of  the  cell  was 
that  the  particle  was  born  in.  Each  particle  represents  electrons  that  were 
born  in  two  cells,  since  they  are  injected  in  alternating  cells.  Therefore, 
we  count  the  current  (and  ionization  rate)  contribution  of  each  particle 
twice.  If  the  current  due  to  a single  particle  is  J,  for  example,  we  assign 
this  value  to  the  cell  in  which  the  particle  stops  at  the  end  of  a time  step. 

In  addition,  the  particle  contributes  an  amount  7 J to  the  current  in  the 

o 

six  adjacent  cells.  Between  this  spreading  and  the  alternation  of  injected 
cells,  a very  smooth  distribution  in  both  space  and  time  can  be  achieved, 
even  when  we  inject  every  other  time  cycle. 

Problems  arise  near  boundaries  because  of  the  lack  of  particles  to 
make  their  1/6  “ contribution.  This  is  partially  compensated  by  allowing  a 
particle  to  penetrate  the  boundary  and  spread  its  charge  back  before  killing 
it.  The  same  is  true  when  a particle  enters  a gamma  thick  pole.  They  are 
not  killed  at  all  when  they  enter  a gamma  thin  pole. 

An  additional  problem  arises  near  the  top  and  front  boundaries, 
where  the  gamma  rays  enter  the  box.  Without  a particle  injecting  boundary 
condition,  we  must  allow  a reasonable  distance  for  the  particle  distribution 
to  build  up.  For  gamma  energies  of  about  1.5  or  2 McV,  a reasonable 
distance  is  1.2  meters.  This  is  illustrated  in  Figure  11.  Here  we  plot  the 
spatial  distribution  of  the  horizontal  current  parallel  to  the  gamma  flux 
(Jz)  as  a function  of  z.  The  time  is  during  the  period  when  the  flux  is 
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z (meters) 


Figure  11.  Comparison  between  current  distribution  calculated  with  parti- 
cles and  ideal  distribution  in  direction  parallel  to  ground 
and  flux.  Time  is  during  rise  of  pulse.  Magnitude  of  "ideal" 
distribution  is  normalized. 


rising  exponentially.  The  dashed  line  indicates  the  manner  in  which  the 
current  should  be  falling  as  a function  of  z (it  is  normalized  and  is  not 
meant  to  show  that  the  magnitude  of  the  current  is  close  to  the  "real" 
magnitude).  Note  the  way  in  which  the  current  builds  up  as  we  move  away 
from  the  front  boundary.  The  flux  angle  of  incidence  is  20°  from  the 
horizontal,  so  there  is  a similar  problem  as  we  move  downward  from  the  top 
boundary;  however,  the  buildup  occurs  over  a distance  of  only  0.4  meters,  in 
this  case,  because  of  the  shallow  incidence  angle. 
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The  fields  predicted  in  this  region  are  complicated  by  the  fact 
that  a radiated  field  generated  by  the  air/ground  interface  is  also  building 


3.4  AIR  CHEMISTRY  EQUATIONS 


From  Section  2.7,  the  air  chemistry  equations  solved  by  DAVID  and 


DAVEJR  are 


IF  “ Se  - (6  - G)ne  - an+ne  , 


IF  = Bne  ' rn+n-  * 


n+  = ne  + n_  , 


ng  = electron  density 

n = negative  ion  density 

n+  = positive  ion  density 

S = ionization  rate 
e 

8 = electron  attachment  rate 
G = electron  avalanching  rate 


a = electron-ion  recombination  rate 
T = ion-ion  recombination  rate. 


These  equations  are  solved  in  a manner  that  is  more  sophisticated 
and  time  consuming  than  necessary,  but  is  consistent  with  the  exponential 
solution  used  in  the  field  and  momentum  equations.  This  will  allow  the  code 
to  use  larger  time  steps,  when  other  factors  allow,  than  could  be  used  with 
simple  differencing. 
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Let  $T  = 8 - G.  Then  eliminating  the  positive  ion  density  from 
Equations  48  and  49.  we  have  the  following  two  equations  to  solve 


U1  e 2 

IF  + (6T  + an-)ne  + ane  = Se  * 


IT  + ^ne)n_  + Fn_  = gne 


These  equations  are  both  of  the  form 


£ + Bf  + Af2  = -D 
dt 


The  solution  of  this  equation  is  given  by  integrating  the  expression 


= - dt  . 


Af  + Bf  + D 


r df 
) Af2  + Bf  + D 


1 „ 2Af  + B -Vq 

X,og  ^ , q > 0 

Vq*  2Af  + B + Vq 


1 _ -1  2Af  + B 

tan  ==— 

V-q 


, q < o 


7 2 

where  q = B - 4AD.  Now  D < 0,  so  that  q > 0.  Using  d = -D,  q = B + 4Ad 


2Af  + B - Vq  _ „ -Vq  t 
2Af  + B + Vq-  " 0 


We  assume  that  the  drivers  (S  and  gn  = d)  are  constant  over  a 

e e 


time  step  At.  Let  f^  be  the  value  of  f just  before  the  new  interaction.  Then 


(jllVslW  ( 

_ \ 2A  /V 


-Vq'  At  (6  - Vq  ' 


1 - K0e"^At 
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where  KQ  = (2AfQ  + B - Vql/(2Afp  + B + VD  • Multiplying  numerator  and 
denominator  by  2A,  yields  the  form  used  in  SUBROUTINE  AIRCHEM: 


f = 


(Vq  + B)K0e'^At  + (Vq  - B) 


2A 


(l- 


K„e 


-Vqit 


) 


(57) 


Table  7 shows  how  these  variables  relate  to  the  air  chemistry  quantities. 
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It  should  be  noted  that  because  ox  numerical  difficulties  in  the 

limit  of  small  electron-ion  densities,  the  density  arrays  must  be  seeded  with 

6 —3  — 3 

an  initial  value  of  10  m (1  cm  ) . A larger  value  may  be  needed  with 

computers  carrying  fewer  significant  figures  than  the  CDC  7600  (60  bit  word) . 

This  point  is  not  very  important,  since  we  initialize  the  densities  the 

first  time  that  they  are  calculated  in  a cell  and  these  initialized  values 

6 “3 

are  generally  in  excess  of  10  m . 


The  initialization  of  the  electron  and  ion  densities  takes  advantage 
of  the  fact  that  the  source  is  rising  exponentially  from  t = It  also 

assumes  that  only  attachment  is  important.  The  initial  values  are 


Table  7.  Relationship  between  general  equation  variables  and 
the  specific  air  chemistry  quantities  (variables  in 
parentheses  are  names  used  in  SUBROUTINE  AIRCHEM). 


Generalized 
Quanti  t,v 

Electron 

Equation 

Negative  Ion 
Equation 

f 

ne(NE(I,J,K)) 

n_(NM(I,J,R)) 

A 

a (ALPHA) 

r(GAM) 

B 

8j  + an_ 

lne 

d 

Se(S(I,J,K)) 

3ne(BETA*NE(I,J,K)) 

60 


j-i -I'i?  j-.  a jji4.  ^ v 'rD~ ^??s>ji3siP5f 


y-^y. 


0 _ e 
e “ B + b * 


n°  = 8:)®/b  , 
— 6 


where  b is  the  e-fold  rate  of  the  rising  Y curve  (t  ~ e ) and  S is  the 

6 

first  value  of  the  ionization  rate  calculated  by  the  particle  routine. 


In  AIRCHEM,  the  electron  equation  is  solved  first  and  the  average 
value  of  the  new  and  old  electron  density  is  used  in  the  driver  function 
for  the  negative  ion  density. 
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SECTION  4 
NUMERICAL  RESULTS 


In  this  section  we  will  review  some  of  the  numerical  calculations 
that  were  used  to  verify  the  operation  of  DAVID  and  DAVEJR.  Many  tests  were 
made  with  the  individual  subroutines  before  they  were  placed  in  the  main  pro- 
gram. The  results  shown  here  are  a sampling  of  those  which  were  made  when  the 
codes  were  fully  assembled. 

When  discussing  the  calculations,  we  will  have  occasion  to  refer 

to  two  sets  of  fields:  transverse  magnetic  (TM)  and  transverse  electric  (TE) . 

The  TM  fields  (E  , E , and  H ) are  those  which  would  exist  in  the  absence  of 
y z x' 

any  object  in  the  mesh.  The  TE  fields  (Ex>  H^,  and  Hz)  are  generated  by  the 
object.  Ideally,  they  are  zero  in  the  absence  of  the  object  and  form  a very 
sensitive  test  of  how  well  boundary  conditions  and  other  aspects  of  the  compu- 
tation are  going.  One  must  be  careful  when  looking  at  parameter  studies  of  TE 
field  behavior  because  they  often  grossly  exaggerate  small  problems. 

4.1  EX  STUDY 

A series  of  investigations  was  made  to  look  at  the  contours  of  Ex 
on  the  air/ground  interface.  Parameters  that  were  varied  are  (1)  the  presence 
and  absence  of  a pole,  (2)  the  angle  of  incidence  (0°  and  20°)  and  (3)  the 
use  of  the  ambient  boundary  condition  compared  to  a perfectly  conducting  outer 
boundary.  The  geometry  and  dimensions  of  the  problem  are  shown  in  Figure  12. 
There  are  seven  0.1  m cells  within  the  mesh  in  the  x-direction  (NX  = 8)  and 
eleven  0.2  m cells  in  the  z-direction  (NZ  = 12).  In  the  air,  the  vertical 
step  size  is  0.2  m (NY  = 9)  while  in  the  ground  the  vertical  step  size  is  0.1  m 
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Figure  12.  Geometry  used  in  EX  study  (only  cells 
near  boundary  are  shown). 

(NYC  = 13).  Ground  drivers  are  used  in  the  first  five  layers  of  ground 
(NRY  = 5).  DAVEJR  was  used  for  these  calculations. 

The  pole  is  four  cells  high  in  the  air  (0.8  m)  and  reaches  down 
to  a perfectly  conducting  plate  at  the  bottom  of  the  ground  mesh  (1.2  m).  It 
is  one  cell  wide  in  the  z-direction  (0.2  m) . Including  its  image  on  the 
opposite  side  of  the  TM  symmetry  plane,  the  pole  is  three  cells  wide  (0.3  m) 


in  the  x-direction.  Remember  that  the  symmetry  plane  passes  halfway  through 
the  first  cell. 

If  the  current  on  the  pole  was  vertical  and  evenly  distributed 
around  its  circumference,  the  Ex  field  that  we  would  expect  to  see  would 
be  the  x-component  of  a radial  electric  field  caused  by  the  accumulation  of 
charge  along  the  pole.  Ex  would  then  peak  in  the  x-direction  and  go  to  zero 
in  the  ± z direction,  behaving  as  cos<j)  in  between,  where  is  measured  from 
the  x-axis. 


If  the  pole  current  flows  around  the  pole  from  front  to  rear  (-z 
to  +z),  we  would  expect  a distribution  of  Ex  which  is  zero  in  both  the  z- 
and  x-directions.  It  would  be  negative  in  one  quadrant  and  positive  in  the 
other  (looking  at  just  one  side  of  the  symmetry  plane).  If  the  pole  has  a 
circular  cross  section,  the  peaks  would  occur  at  4>  = ± 45°,  but  with  a 
rectangular  cross  section,  the  distribution  will  be  shifted. 

The  figures  that  follow,  show  contours  of  constant  Ex  computed  at 
the  air-ground  interface.  They  were  hand  drawn  from  printed  data  and  so  are 
not  very  accurate.  The  general  distributions  should  be  representative,  however. 
All  of  the  figures  have  one  error,  which  does  not  affect  our  analysis.  The 
contours  are  drawn  in  such  a way  that  the  field  goes  to  zero  at  the  side  of 
the  pole.  Ex  does  not  do  this,  of  course,  but  should  appear  to  eminate  from 
here  (+x  side) . It  does  go  to  zero  on  the  front  and  back  (-z  and  +z  sides) . 

Figures  13a  - 13d  show  the  Ex  contours  which  are  generated  when  the 

pole  is  excited  by  a source  incident  horizontally  (9  = 0°).  The  peak  dose 
13  -8 

rate  is  10  rad/sec  and  peaks  near  4 x 10  second.  The  contours  are  shown 

at  5 x 10~9,  5.8  x io"9,  9.8  x 10"9,  and  2.9  x io"8  second.  The  ground 
_2 

conductivity  is  10  “ mho/m.  The  "ambient"  boundary  condition  is  used  (see 
Section  3.1).  The  field  contours  are  clearly  characteristic  of  a current 
running  around  the  side  of  the  pole  at  early  times  (Figures  13a, b).  Note  that 
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the  contours  are  labeled  in  units  of  volts/meter.  The  numbers  in  the  upper 
corners,  i.e.,  JR22.17(1)  in  this  case,  denote  the  computer  run  number  so 
that  we  can  refer  to  it  in  the  future.  At  later  times  (Figure  13c, d),  the 
negative  fields  generated  near  the  outer  x boundary  begin  to  influence  the 
fields  near  the  pole.  The  influence  is  not  great  enough  to  seriously  alter 
the  calculation  of  currents  on  the  pole.  Another  effect  becomes  perceptible. 
The  distribution  shifts  to  the  left  indicating  the  influence  of  a net  charge 
accumulation  near  the  ground. 

In  Figure  14,  we  show  the  fields  generated  when  the  gammas  are 
incident  at  an  angle  of  20°.  Figures  14  a and  14b  are  at  the  same  times  as 
Figures  13c  and  13d.  In  this  case,  the  presence  of  a charge  on  the  pole, 
caused  by  a net  vertical  flow  of  current,  is  quite  obvious.  The  positive 
field  contours  are  barely  influenced  by  the  outer  boundary  and  are  shifted 
"downstream"  by  the  added  contribution  due  to  the  current  flowing  around  the 
circumference.  The  "downstream  shift"  is  even  more  obvious  when  the  pole 
is  gamma  thick  (opaque  to  gamma  rays),  as  is  seen  in  Figure  i5. 
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The  effect  of  using  perfectly  conducting  walls,  instead  of  the 

"ambient"  boundary  condition  is  shown  in  Figure  16.  The  positive  fields 

generated  by  the  pole  are  now  confined  to  a small  region  near  the  pole.  The 

influence  on  the  pole  currents  is  not  as  strong  as  might  be  indicated  by  the 

contours,  but  is  still  significant  with  the  boundaries  this  close.  The 

variation  of  H in  the  z-direction,  under  the  influence  of  both  the  ambient 
x 

and  conducting  boundary  conditions  is  shown  in  Figure  17. 

Figure  18  shows  the  Ex  contours  which  are  generated  by  DAVEJR  in 
the  absence  of  a pole,  for  comparison  with  the  previous  figures.  In 
the  region  of  the  pole  location,  the  fields  are  generally  negligible.  Ex 
should  not  exist  at  all.  Closer  inspection  of  the  computer  data  shows  that  it 
is  generated  mainly  by  the  poor  outer  boundary  condition  in  the  ground;  it 
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Figure  14b.  Approximate  contours  of  constant  Ex  at  T = 2.90E-8  sec 
(Ambient  boundary  condition,  0 = 20°,  Y = YAMAX). 


decreases  as  we  look  higher  in  the  air.  In  fact,  very  large  TE  fields  are 
generated  deep  in.  the  ground.  These  are  still  small  compared  to  the  TM 
fields,  but  influence  them  in  a non-negligible  manner.  More  work  should 
be  done  on  improving  the  ambient  boundary  condition  in  the  ground. 


Figure  17. 
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Variation  of  HX  in  the  vicinity  of  the  pole,  comparing  the 
field  seen  with  the  ambient  boundary  condition  and  the 
perfectly  conducting  boundary  condition.  The  variation  is 
along  a line  0.1m  above  the  ground  and  passing  through  the 
pole  in  the  z-direction. 
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Figure  18a.  Approximate  contours  of  constant  Ex  at  T = 9.3r-9  sec. 

(Ambient  boundary  condition,  0 = 20°,  Y = YMAX,  no  pole) 
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Figure  18b.  Approximate  contours  of  constant  Ex  at  T = 2.9E-8  sec. 

(Ambient  boundary  condition,  0 = 20°,  Y = YMAX,  no  pole). 
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4.2  GLANC/DAVID  COMPARISON 


In  order  to  make  sure  that  DAVID  was  behaving  in  a reasonable  manner, 
we  made  a comparison  between  DAVID,  without  an  object,  and  the  one-dimensional 
close-in  coupling  code,  GLANC2.  There  are  several  ways  in  which  we  could  have 
performed  this  experiment.  We  chose  to  make  the  GLANC  calculation  as  close  as 
possible  to  the  DAVID  calculation,  i.e.,  the  same  vertical  step  size,  etc. 

In  this  way  we  can  see  if  DAVID  gives  as  reasonable  an  answer  as  GLANC  does 
under  the  same  circumstances;  however,  because  the  grid  sizes  are  larger  than 
would  normally  be  considered  practical  in  a pure  environment  calculation,  we 
will  not  be  comparing  the  DAVID  prediction  to  the  "real"  environment. 


The  question  of  spatial  resolution  is  one  which  must  be  addressed  at 
some  point.  Just  exactly  how  much  resolution  is  required  to  do  the  calculation 
properly?  The  DAVID  geometry  used  in  this  comparison  is  the  same  as  in  the 
previous  (EX)  study.  Therefore,  each  cell  is  20  cm  high  (DY)*,  20  cm  long  (DZ) , 
and  10  cm  wide  (DX) . Now,  there  is  a significant  amount  of  field  variation 
occurring  within  20  cm  of  the  ground  in  a close-in  environment  calculation. 

The  same  would  be  true  near  the  surface  of  an  object.  It  is  important  to  know 
whether  we  are  simply  losing  space  and  time  resolution  or  whether  the  large 
cell  sizes  are  causing  us  to  calculate  a totally  incorrect  answer.  Our  studies 
indicate  that  we  are  simply  losing  resolution.  There  have  not  been  enough 
tests  to  fully  confirm  that  fact  however.  Assuming  that  our  problem  is  basically 
one  of  resolution,  what  is  the  impact? 

DAVID  is  intended  to  calculate  voltages  and  currents  on  objects 
within  the  deposition  region.  The  object  is  going  to  be  resolved  to  the  same 
degree  as  the  fields  that  we  are  coupling  to  it.  The  fact  that  we  may  not  be 
resolving  the  short  wavelength  components  of  the  fields  is,  at  least,  consistent 


* 10  cm  in  the  ground. 


■*  ! 


with  the  fact  that  we  could  not  calculate  the  response  of  the  object  to 
those  wavelengths  anyway.  Therefore,  there  seems  to  be  little  reason  to 
worry  about  them. 

A few  words  about  1-D  calculations  are  also  appropriate.  Consider 
a plane  wave  of  gamma  rays  sweeping  across  'the  surface  of  the  earth  in  the  y-z 
direction  with  an  angle  of  incidence  0 (see  Figure  6) . In  the  absence  of  an 
object,  there  is  no  variation  of  flux  or  fields  in  the  x-direction.  Also, 
if  the  attenuation  of  the  gamma  rays  can  be  ignored  in  the  z-direc*’?r.,  over 
a distance  which  is  important  to  the  solution  of  the  problem  (limited  by 
air  conductivity),  the  physics  observed  at  any  point  in  the  z-direction  will 
be  the  same  as  at  any  other  point,  except  that  the  time  history  will  be 
delayed  by  the  amount  of  time  that  it  takes  light  to  reach  the  observer  from 
the  burst.  Therefore,  derivatives  with  respect  to  z are  related  to  the  time 
derivative  through  the  speed  of  light  and  the  cosine  of  the  angle  of  incidence. 
In  the  absence  of  a ground  plane,  the  same  would  be  true  of  a derivative  with 
respect  to  y,  except  that  it  is  proportional  to  sin0.  The  introduction  of  a 
ground  plane  introduces  additional  contributions  to  the  y-derivative. 

Under  these  circumstances,  one  can  calculate  the  fields  using  only 
the  y-derivatives,  i.e.,  the  problem  can  be  solved  by  differencing  along  a 
vertical  line.  This  is  a realistic  thing  to  do  when  calculating  fields  in 
the  close-in  region  at  times  when  the  conductivity  is  high  enough  to  limit 
the  distance  over  which  an  observer  can  see  to  those  which  are  small  compared 
to  a gamma  mean  free  path  (~100  m) . The  GLANC  code  is  based  on  this 
principle. 

We  will  now  briefly  derive  a set  of  1-D  equations  for  the  purpose 
of  (1)  illuminating  the  principles  involved  and  (2)  deriving  a set  of  equa- 
tions which  will  help  to  interpret  and  check  the  numerical  calculations  made 
with  GLANC  and  DAVID.  The  equations  we  will  derive  are  not  in  the  same 
form  as  used  by  GLANC,  but  that  will  be  of  no  consequence. 


We  start  with  the  field  equations  in  our  units  (h  = ZqH,  3 = Z^J, 


Zq  = 120tt  ohms,  see  Section  2.2) 
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(60a) 


(60b) 


Only  the  TM  mode  will  be  excited  by  the  currents  which  are  in  the  y-z 
direction.  Therefore,  only  the  hx,  E^,  and  Ez  fields  need  be  considered. 
Remembering  that  derivatives  with  respect  to  x are  zero,  we  have  three  scalar 
equations 
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The  z-derivative  is  simply 
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(62a) 


The  y-derivative  has  an  extra  term,  since  variations  in  the  y-direction  are 
caused  by  both  the  time  phasing  and  a discontinuity  of  the  medium,  i.e., 
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The  primed  variable  (y*)  reminds  us  that  we  have  removed  the  time  phasing 
effect. 
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Equations  61  become,  after  rearrangement, 
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These  are  equivalent  to  those  used  in  GLANC.  Under  conditions  of  high 
conductivity  when  displacement  current  is  much  less  than  the  conduction 


current  (3E/3t  « — E),  the  first  two  can  be  written 
e0 
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Note  that  we  have  returned  to  the  true  current  density  J rather  than  the 


normalized  current  density  3"  = 


Equation  63d  will  be  useful  in  checking  the  numerical  calculations*. 
It  states  that  the  vertical  E-field  is  generated  by  both  the  vertical  current 
and  the  time  derivative  of  the  magnetic  field.  There  is  no  explicit  spatial 
derivative.  The  horizontal  field  equation  does  contain  a y' -derivative  and 
is  therefore  more  difficult  to  use  in  checking  numerical  results. 


It  will  also  be  useful  to  understand  the  relation  between  the 
electric  and  magnetic  fields  at  early  times,  when  the  conduction  current  is 
negligible  compared  to  the  displacement  current.  Under  these  circumstances. 
Equations  63a, b reduce  to 


* A peak  Y of  10^  rad/sec  is  used,  as  in  the  Ex  study.  Therefore,  much 
of  the  calculation  will  be  for  times  when  this  approximation  is  valid. 
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The  GLANC/DAVID  comparisons  are  shown  in  Figures  19  through  23. 

It  was  found  to  be  more  practical  to  make  side-by-side  comparisons  rather 
than  overlays.  Since  t = 0 is  defined  differently  in  the  two  codes,  we 
have  drawn  a dashed  line  in  each  graph  which  shows  the  time  that  the  t time 
history  peaks.  Note  that  the  GLANC  calculations  extend  20  nsec  past  the  y 
peak,  while  the  DAVID  calculations  extend  only  5 nsec. 

The  fields  represent  those  which  would  be  seen  10  cm  above  the 

ground.  H , E , and  the  sources  are  actually  calculated  at  that  height, 
x y 

GLANC  interpolates  E z to  get  it  there,  while  in  DAVID  we  simply  output  the 

quantities  corresponding  to  the  first  cell  above  the  ground.  E is  therefore 

z 

the  value  calculated  at  20  cm,  as  is  E . H is  calculated  at  20  cm,  also, 

x y 

while  t!,  is  computed  at  10  cm.  The  latter  three  field  components  should 
ideally  be  zero.  They  are  not,  of  course,  and  their  amplitude  relative  to 
the  TM  mode  fields  is  an  indication  of  when  the  calculation  is  starting  to 
deteriorate.  We  reiterate  the  fact  that  the  presence  of  an  object  within 
the  grid  will  greatly  control  the  behavior  of  the  TE  mode  fields.  Therefore, 
the  fact  that  relatively  large  spurious  TE  fields  develop  in  this  calculation, 
after  the  gamma  peak  does  not  mean  that  the  problem  will  develop  on  the  same 
time  scale  when  an  object  is  present. 

GLANC  and  DAVID  predict  about  the  same  ionization  rate  at  the  time 
of  the  y peak  (Figure  19),  but  this  is  not  representative  of  the  overall 
agreement.  There  is  more  structure  in  the  DAVID  curve  so  that  it  is  some- 
times greater  than  the  GLANC  prediction  and  sometimes  less.  The  peak  ioniza- 
tion rate  predicted  by  DAVID  is  1.5  times  greater  than  GLANC  while  at  10  nsec 


before  the  '■  peak,  it  is  1.5  times  less.  Both  codec  injected  particles 
every  time  step  for  this  comparison.  Normally,  we  inject  DAVID  every 
other  tir.c  step  in  order  to  reduce  the  number  of  particle".  At  very  early 
times  (about  two  orders  of  magnitude  down  from  the  peak)  the  sources  calculated 
with  injection  every  other  time  step  can  be  as  much  as  40  percent  off  those 
calculated  with  injection  every  time  step.  However,  as  the  number  of  particles 
builds  up  and  the  rise  rate  of  the  t curve  decreases,  the  two  sources  agree 
much  better.  We  ran  this  same  DAVID  calculation  with  alternate  cycle  injec- 
tion and  found  5 percent  agreement  at  the  source  peak. 

The  onductivity  (Figure  20)  predicted  by  DAVID  at  the  ^ peak  is 
1.5  times  greater  than  that  seen  in  GLANC.  Ten  nanoseconds  before  this  it  is 
1.7  times  lower.  The  DAVID  calculation  did  not  go  late  enough  in  time  to 
compare  the  peak  conductivities. 

DAVID  will  give  lower  source  strengths  during  the  y rise  for  ob- 
servers in  the  first  cell  above  the  ground  because  the  spatial  averaging 
scheme  (Section  3.3)  does  not  account  for  the  fact  that  the  cell  in  the 
ground  does  not  contribute  its  share  to  the  sources  in  the  cell  just  above 
the  ground.  This  can  easily  be  fixed,  but  was  not  considered  to  be  an 
important  enough  problem  to  warrant  the  more  complicated  logic  that  would 
be  required  in  the  initial  version  of  the  code. 

In  the  Compton  current  comparisons  (Figure  21)  we  start  to  see 
some  interesting  differences,  which  are  fortunately  explainable.  The  horizontal 
currents  agree  to  the  same  degree  as  the  ionization  rate,  but  GLANC  shows  a 
vertical  current  with  far  more  fine  structure  5 to  10  nsec  before  the  t 
peak.  The  GLANC  current  actually  changes  sign  for  a short  time.  This  is 
not  a physically  real  effect  and  would  not  ordinarily  be  predicted  by  the 
code.  It  is  caused  by  the  fact  that  we  forced  GLANC  to  inject  particles 
parallel  to  the  gamma  path.  Since  they  are  all  traveling  in  the  same  direc- 
tion, the  effect  of  magnetic  turning  is  overemphasized.  Normally,  GLANC 
would  inject  a distribution  of  particles  at  angles  about  the  ray.  The  center 
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forces  as  the  individual  particles  themselves.  In  DAVID  the  spatial  averag- 
ing also  acts  as  a time  averaging  scheme.  Half  of  the  current  scored  in  a 
given  cell  is  contributed  by  the  particles  within  the  cell.  The  other  half 
is  contributed  by  the  particles  in  the  six  adjacent  cells  (the  contribution 
from  one  cell  is  missing  at  the  boundaries).  Since  the  particles  in  the 
adjacent  cells  are  at  different  phases,  there  is  a time  averaging  effect. 

Note  that  J is  two  orders  of  magnitude  lower  than  J and  J , so  it  is  not 
x y z 

a significant  perturbation. 

The  horizontal  electric  fields  (Ez)  predicted  by  the  two  codes 
(Figure  22)  agree  well,  within  a factor  of  1.5.  The  early  time  vertical 
fields  do  not.  In  fact,  they  have  the  opposite  sign  until  5 nsec  before  the 
Y peak.  The  DAVID  fields  maintain  the  same  sign  throughout  the  time  history. 
GLANC  starts  out  with  the  opposite  sign  and  then  changes  it  to  agree  with 
DAVID.  The  peak  vertical  fields  agree  to  within  a factor  of  1.5 

The  early  time  discrepancy  in  the  vertical  fields  is  undoubtedly 
due  to  the  way  in  which  the  propagated  field  is  handled  in  each  code.  Neither 
code  does  very  well,  actually,  during  this  low  conductivity  phase  of  the 
problem.  DAVID  tends  to  underestimate  the  propagated  part  of  the  signal  be- 
cause the  sources  and  fields  are  built  up  over  a relatively  short  distance 
fr<vn  the  front  boundary.  GLANC  assumes  that  the  sources  are  the  same,  as 
a f met ion  of  local  time  all  the  way  back  to  z = This  would  tend  to 
givt  a larger  propagated  signal  than  DAVID.  The  1-D  approximation  is  not 
as  tad  here  as  one  jnight  first  expect.  The  fact  that  the  sources  are  assumed 


to  bo  the  same  back  to  z = •»  is  really  not  important  when  the  burst  is  off 
the  ground  because  the  gamma  wave  front  arrives  faster  than  the  electric  field, 
which  is  generated  along  the  air/ground  interface.  Therefore,  conductivity  is 
allowed  to  build  up  and  absorb  the  propagated  field.  The  field  would  not  grow 
large  in  any  case,  because  it  loses  phase  with  the  Compton  current.  Also, 
the  longer  integration  distance  partially  compensates  for  the  fact  that  the 


sources  are  stronger  near  ground  zero'. 


The  1-D  approximation  becomes  worse  as 
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the  angle  of  incidence  approaches  zero.  In  our  problem  the  GLANC  prediction 
is  probably  reasonable.  We  will  return  to  this  point  after  discussing  the 
magnetic  field  predictions. 

The  magnetic  fields  are  shown  in  Figure  23.  GLANC  outputs  its 
magnetic  fields  in  units  of  gauss,  so  an  extra  volts/meter  scale  has  been 
added  to  facilitate  comparison  with  DAVID*.  The  GLANC  predictions  are  a 
factor  of  three  larger  than  the  DAVID  predictions.  An  inspection  of  the 
spatial  plots  revealed  that  the  problem  was  caused  by  not  having  the  grid 
large  enough.  One  could  see  the  H-field  building  up  in  the  +z  direction, 
even  on  the  rise  of  the  pulse  when  the  field  should  be  decreasing  as  a 
function  of  increasing  z. 

This  brings  out  a very  important  point,  which  should  be  remembered 
when  performing  3-D  calculations.  Magnetic  fields  are  not  nearly  as  intimi- 
dated by  conductivity  as  electric  fields.  The  skin  depth  arguments  which 
are  commonly  used  to  claim  that  the  conductivity  is  high  enough  to  isolate  an 
object  from  its  surroundings  do  not  apply  so  well  when  magnetic  field  coupling 
is  involved.  This  is  emphasized  by  the  later  time  DAVID  calculations  which 
show  H and  H growing  substantially  and  H falling  rapidly  because  of  this 

/ Z X 

boundary  problem.  The  conductivity  at  this  time  is  very  high  (—0.5  mho/m). 

We  considered  this  problem  before  doing  the  pole  calculations  in 
the  next  section.  The  grid  was  enlarged  in  the  z-  and  y-directions.  It  isn't 
that  critical,  however,  as  the  spatial  plots  will  show.  After  the 
conductivity  builds  up  and  starts  to  control  the  current  running  on  the 
vertical  pole,  the  magnetic  field  generated  by  the  ground  is  seen  to  be 
dominated  by  that  caused  by  the  current  running  on  the  pole.  The  pole  coup- 
ling in  this  case  is  more  electric  than  magnetic. 


* H(volts/meter)  = 3 x 104B(gauss). 
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We  now  return  the  effect  of  the  propagated  field  on  the  early  time 
E prediction.  We  know  that  this  (E^)  field  is  composed  of  two  parts:  the 

vertical  component  of  the  local  radial  field  and  the  part  that  has  propagated 
along  the  earth.  In  order  to  give  the  Poynting  vector  (E  x il)  the  right 
direction  (+z) , the  propagated  Ey  component  must  be  negative  and  this  will 
then  oppose  the  locally  generated  component.  This  is  shown  by  Equations 
63f  and  63g,  which  are  zero  conductivity  approximations. 

If  there  were  no  propagated  fields,  both  E and  E would  be  positive 

y z 

in  our  coordinate  system.  The  magnetic  field,  which  is  positive  during  these 
times,  would  therefore  decrease  Ey  and  increase  Ez  (remember  that  if  there 
was  no  discontinuity  to  generate  Hx>  there  would  also  be  no  propagated  signal). 

A little  work  shows  that  both  the  DAVID  and  GLANC  calculations  satisfy  these 
equations  at  early  times  (the  y' -derivative  term  in  E can  be  ignored  as  a first 
approximation) . The  propagated  field  contribution  is  simply  much  larger  in  GLANC, 
causing  a sign  reversal  at  early  times,  until  the  conductivity  is  large  enough 
to  make  Equation  65d  valid.  At  times  near  and  after  the  peak,  the  electric 
fields  are  behaving  as  J/o  in  both  calculations  and  therefore  agree  much  better. 

4.3  POLE  STUDY 

Having  gained  confidence  in  DAVID  through  the  GLANC  comparison 

descij...^d  in  the  previous  section,  and  feeling  that  we  understood  most  of 

the  code's  problems  and  limitations,  calculations  were  performed  for  the 

problem  of  a vertical  pole  passing  through  the  earth  and  exposed  to  a peak 
13 

Y of  10  rad/sec.  An  angle  of  incidence  of  20°  was  chosen  because  a few 
experimental  data  points  were  available  to  tell  us  whether  the  current 
induced  on  the  pole  was  a reasonable  prediction.  The  problem  geometry  is 
shown  in  Figure  24. 

Several  parameters  were  varied,  including  the  width  of  the  pole  in 
the  x-direction.  The  data  had  been  taken  for  a pole  with  a circular  cross 
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Seometry  used  in  the  pole  current  calculation.  Note  the  use  of 
two  pole  sizes.  The  fat  pole  has  a width  (Ax)  of  30  cm  and  a 
circumference  of  100  cm.  The  narrow  pole  has  a width  of  10  cm 
and  a circumference  of  60  cm.  The  dimensions  include  the  image 
half  of  the  pole  on  the  other  side  of  the  symmetry  plane.  Cells 
are  shown  only  near  the  boundaries. 


section  and  a diameter  of  10  cm.  With  our  constant  spatial  step  sizes,  we 
could  not  calculate  a pole  that  narrow  and  still  maintain  a grid  large 
enough  to  keep  the  outer  boundary  a sufficient  distance  away.  We  felt, 
however,  that  during  the  highly  conducting  portion  of  the  signal,  that  we 
could  scale  the  pole's  surface  currents  by  its  circumference.  Therefore 
we  ran  parameter  variations  using  two  pole  sizes:  100  cm  circumference 

and  60  cm  circumference.  The  10  cm  circular  pole  has  a circumference  of 
31.4  cm.  After  scaling  the  net  currents  by  the  ratio  of  the  pole  circum- 
ference to  31.4  cm,  we  found  that  the  answers  agreed  quite  well  in  the  air. 
The  pole  current  at  a depth  of  -45  cm  was  found  to  be  fairly  insensitive  to 
pole  size.  These  studies  will  be  shown  shortly. 


The  scaling  in  the  air  was  expected  to  work  because  the  poles  were 
fairly  close  in  size  and  their  cross  sectional  shapes  were  not  extremely 
different.  The  near  fields  of  the  pole  should  vary  logarithmically  with 
effective  pole  radius.  As  we  shall  see,  the  currents  running  on  each  face 
of  the  pole  can  vary  significantly.  It  would  have  been  better  to  scale  each 
component  separately  before  summing,  instead  of  scaling  the  net  current. 
However,  the  latter  method  worked  well  enough  for  our  purposes. 


The  reasons  for  the  insensitivity  of  the  pole  currents  below 

ground  to  pole  size  will  be  discussed  later.  Most  of  the  calculations  were 

performed  using  ground  drivers*,  which  fell  off  exponentially  with  the 

slant  range  into  the  ground.  Ground  conductivity  enhancement  was  sometimes 

used.  The  enhanced  conductivity  was  assumed  to  be  proportional  to  the 

-14 

local  gamma  flux,  with  the  constant  of  proportionality  being  10 

13 

mho/m/ rad/sec . Thus,  with  a peak  gamma  flux  of  10  rad/sec,  the  maximum 

enhancement  at  the  surface  would  be  0.1  mho/m.  The  ambient  ground  conductivity 
-7  _ 

was  set  at  10  ' mho/m.  lhe  dielectric  constant  was  10. 


* When  viewing  the  curves,  it  can  be  assumed  that  ground  currents  are 
present  unless  otherwise  stated. 


Three  steps  are  required  to  calculate  the  pole  current:  (1)  the 

tangential  magnetic  fields  are  extrapolated  to  the  pole  surface,  (2)  the 
surface  current  density  at  that  point  is  computed,  and  (3)  the  current  density 
is  integrated  around  the  pole.  Since  DAVID  is  a general  purpose  code, 
intended  to  handle  many  • '•'tries,  it  was  felt  that  it  would  be  most 

efficient  to  use  a separat  .urrent  calculation  algorithm  for  each  type  of 
object  rather  than  build  in  a general  one,  which  would  require  a great  deal 
of  computer  logic.  Since  the  code  requires  the  user  to  supply  his  own 
algorithm,  we  will  go  into  more  detail  on  how  we  did  it  for  the  pole  than 
would  normally  be  considered  prudent.  There  are  several  peculiarities 
of  the  code  that  may  not  be  obvious  at  first. 

A two  step  extrapolation  is  required  because  the  tangential  magnetic 
fields  are  computed  on  the  cell  walls  and  not  at  the  cell  center  (see  Figure 
8).  In  Figure  25,  we  show  a cross  section  of  the  "fat"  pole  and  the  sur- 
rounding cells  which  are  needed  for  the  extrapolation.  The  pole  is  located 
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Figure  25.  Geometry  used  to  compute  current  running  on  the 
"fat"  pole  (3  cells  wide  in  x-direction, 
including  image). 


at  K = 9,  i.e.,  the  pole  is  centered  at  z = 1.7  m from  the  front  boundary. 

In  the  x-direction,  the  pole  extends  over  I = 1,2.  However,  the  symmetry 
plane  runs  through  the  center  of  the  1=1  cells,  so  half  the  width  of 
that  cell  is  associated  with  the  image  half  of  the  pole.  The  locations  of 
the  fields  that  we  need  for  ,he  pole  current  calculation  are  shown  by 
arrows.  The  first  step  is  to  interpolate  in  order  to  estimate  the  values 
of  the  field  at  the  cell  center  (designated  by  the  open  circles  in  the 
figure).  At  the  symmetry  plane,  we  have  HX(1,J,K)  = HX(2,J,K),  so  no 
interpolation  is  required  as  long  as  we  only  use  linear  interpolations.  Then, 
an  extrapolation  is  required  to  estimate  the  tangential  field  at  the  pole 
surface.  For  these  calculations,  we  used  a linear  extrapolation*. 

Having  obtained  values  of  HX  and  HZ  at  ;.he  centers  of  the  cell 
faces  which  define  the  pole,  the  components  must  be  summed  in  the  proper 
sense  to  calculate  the  net  current.  At  some  point,  the  magnetic  fields 
calculated  by  the  code  must  be  divided  by  ZQ  = 120tr  to  obtain  units  of  amp/m. 
Then,  the  net  current  is  given  by 

I = £ H • d£ 

We  define  positive  current  as  the  flow  of  positive  charge  in  the  +y-direction 
(electron  flow  upward).  Therefore,  the  net  current  is 

I = 5Z  AX-HX  - £ AZ-HZ  - £ AX-HX  + J2  AZ-HZ  (65) 

BACK  SIDE1  FRONT  SIDE2 

where  these  are  the  sums  of  the  renormalized  H-fields  on  each  side.  By 
symmetry,  • 

HZ  lsiDE2  = ' HZlsiDEl  ' 


* Extrapolating  the  product  of  H and  the  cylindrical  radius  from  the  pole 
works  better  in  this  case,  because  of  the  pole  shape.  We  chose  to 
extrapolate  the  H-field  along  because  of  the  more  general  applicability 
of  the  procedure  and  because  it  seemed  to  work  well  enough. 
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Figure  26  shows  the  currents  calculated  on  the  "fat"  pole  by 
DAVEJR  (prescribed  current  and  ionization  rate).  We  also  show  the 
conductivity  time  history  near  the  pole  for  the  30  cm  observer.  No  enhanced 
conductivity  is  used  in  this  calculation.  The  numbers  in  parentheses  near 
the  waveform  peaks  show  the  value  of  the  current  (in  kiloamps)  normalized 
by  the  circumference  of  the  10  cm  diameter  circular  pole.  We  will  state 
the  measured  peak  values  after  completing  the  theoretical  parameter  study. 

The  pole  useu  in  this  and  the  following  calculations  is  transparent  to  gamma 
rays. 

The  first  notable  feature  of  Figure  26  is  the  way  in  which  the 
currents  at  +30  cm  and  +70  cm  follow  the  conductivity.  Also,  they  are 
very  close  in  magnitude  even  though  one  observer  is  near  the  ground  and  the 
other  is  near  the  end  of  the  pole.  Actually,  the  +30  cm  observer  is  not 
that  close  to  the  ground  in  as  much  as  it  is  isolated  by  the  air  conductivity. 
However,  this  isolation  is  not  as  great  as  one  might  guess  at  first.  What 
is  happening  is  that  the  boundary  condition  provided  by  the  pole  is  so  strong 
that  it  controls  the  fields  in  the  poles  vicinity,  even  near  the  ground. 

One  would  probably  get  similar  answers  for  the  current  running  on  the  pole 
in  the  air  without  including  the  ground  plane  at  all*.  Later  studies  will 
show  that  this  is  not  true  for  grazing  incidence  angles  because  of  the 
importance  of  the  propagated  field  generated  at  the  surface.  Note  the 
slight  negative  excursion  at  early  times.  This  will  be  seen  to  be  a 
scattered  field  response.  It  is  not  very  important  in  this  case,  being  a 
factor  of  430  below  the  peak,  but  it  is  of  theoretical  interest.  We  will 
watch  its  behavior  later  as  we  vary  the  flux  and  angle  of  incidence.  This 
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* It  would  not  be  surprising  to  find  that  with  this  incidence  angle  (or  { 

greater)  one  can  obtain  reasonable  answers  using  an  infinite  cylinder  ? 

in  air.  j 
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Figure  26.  Currents  generated  on  fat  pole  without  ground  conductivity  en' 
hancement.  Numbers  in  parentheses  are  peak  scaled  currents 
(kamps).  Air  conductivity  is  also  shown.  Analytic  sources. 
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part  of  the  response  is  more  significant  for  observers  below  the  ground. 

The  current  at  +70  cm  leads  that  at  +30  cm  at  early  times  by  the 
difference  in  time  required  by  the  gamma  wave  front  to  reach  the  two  observers 
(~  1 nsec  with  this  20°  incidence  angle).  The  vertical  arrows  at  the  top  of 
the  figure  show  the  times  at  which  the  t pulse  peaks  at  the  two  observers. 

The  current  at  +70  cm  is  limited  to  a value  less  than  that  at  +30  cm, 
probably  because  there  is  some  charge  build-up  at  the  top,  despite  the  high 
conductivity. 

The  current  seen  at  -45  cm  has  a much  different  character  than 
that  seen  at  positions  in  the  air.  It  will  also  prove  to  be  much  more 
sensitive  to  various  physical  factors,  e.g.,  ground  drivers  and  conductivity 
enhancement.  The  current  starts  out  with  a negative  swing,  which  is  the 
scattered  field  response  (the  electric  field  in  the  ground  is  negative  during 
almost  the  entire  time  frame  plotted).  It  then  swings  positive,  rises  to  a 
peak  and  oscillates.  The  oscillation  is  riding  on  a rising  base.  It  is 
possible  that  this  late  time  rising  is  a numerical  problem  caused  by  the 
outer  boundary  in  the  ground.  However,  it  is  not  unreasonable  to  expect  the 
current  to  try  and  make  a more  uniform  distribution  over  the  pole  at  late 
times,  matching  the  current  on  the  pole  in  the  air,  and  that  is  what  it  ap- 
pears to  be  doing.  We  will  look  at  how  different  physical  effects  influence 
the  current  at  -45  cm  later. 

Figure  27  shows  the  same  calculations  for  the  "narrow"  pole  (60 
cm  circumference  instead  of  100  cm).  Qualitatively,  they  are  the  same. 

When  the  air  observer  currents  are  scaled,  the  peak  currents  are  seen  to  be 
nearly  the  same  as  the  scaled  fat  pole  currents.  That,  of  course,  is  the 
point  we  are  trying  to  prove.  The  same  is  not  true  at  the  ground  observer. 

As  a matter  of  fact,  the  positive  current  becomes  larger  rather  than  smaller. 
The  negative  swing  is  smaller,  as  one  might  expect  from  a scattered  field 
response.  The  negative  peak  is  probably  controlled  more  by  the  rise  of  the 
positive  component  than  by  the  inductance  of  the  pole,  however.  The 
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positive  current  corresponds  to  electrons  running  up  the  pole.  This  current 
is  driven  by  the  voltages  which  are  trying  to  pull  electrons  back  out  of  the 
ground.  The  pole  offers  itself  as  a convenient  path  and,  in  addition,  "amplifies" 
the  voltage.  A narrower  pole  causes  a greater  radial  field  concentration  in 
the  air  and  a larger  current  to  the  surface.  The  current  leaves  the  pole  just 
above  the  surface  because  of  the  hign  air  conductivity  and  the  currents  on 
the  pole  higher  up  do  not  see  any  difference. 

The  effect  of  introducing  ground  conductivity  enhancement  is 
shown  in  Figure  28.  The  currents  running  on  the  pole  in  the  air  are  not 
changed  and  are  therefore  not  plotted.  Only  the  currents  seen  at  -45  cm  on 
the  narrow  pole  are  shown.  The  scattered  field  response  is  not  altered  be- 
cause the  enhanced  conductivity  has  not  risen  high  enough  yet.  It  is  not 
clear  that  it  would  be  anyway,  since  the  layer  of  significant  enhancement 
lies  above  a depth  of  20  cm  and  the  fields  that  excite  the  pole  at  this  time 
can  propagate  down  from  that  layer  faster  than  the  conductivity  could 
change  +hem.  The  main  effect  of  the  enhancement  is  to  dampen  the  positive 
part  of  the  current  time  history.  The  positive  peak  is  lowered  and  the 
oscillations  die  away  faster. 


Until  now,  we  have  only  looked  at  prescribed  current  calculations 
(DAVEJR).  Figure  29  shows  that  the  particle  calculations  (DAVID)  are  not 
significantly  different.  Here  we  compare  the  current  predicted  at  +30  cm 
and  -45  cm.  The  biggest  difference  is  in  the  air,  where  the  currents  build 
up  faster  with  the  prescribed  currents.  The  prescribed  current  formulation 
assumes  that  the  electron  current  is  in  equilibrium  with  gamma  flux.  In 
sea  level  air,  the  lifetime  of  the  Compton  electron  is  on  the  order  of  10 
second.  This  is  at  least  comparable  to  the  time  over  which  the  gamma  pulse 
changes  significantly.  In  this  case,  it  can  be  considered  much  greater 
than  that  time.  Therefore,  an  equilibrium  condition  does  not  exist  and 
electrons  born  during  an  early  part  of  the  gamma  pulse  exist  at  the  same 
time  as  those  born  at  a later  time.  During  the  rise  of  the  gamma  pulse. 
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when  it  is  behaving  exponentially,  this  deviation  from  equilibrium  manifests 
itself  as  a simple  delay  in  the  current.  After  the  peak,  when  the  flux 
varies  more  slowly,  the  two  predictions  agree  better.  If  we  had  put  a 
proper  delay  into  the  prescribed  sources,  the  comparison  would  have  been 
better. 


i 

i 
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The  effect  of  the  delay  is  also  seen  in  the  early  time  current 
induced  on  the  pole  at  -45  cm.  The  current  at  this  time  is  driven  by  the 
vertical  electric  field,  which  is  created  by  the  gradient  of  the  magnetic 
field,  which  has  diffused  downward  from  the  surface.  No  such  delay  is  seen 
in  the  positive  part,  of  the  current.  This  portion  is  highly  influenced 
by  the  ground  drivers  which  are  prescribed  in  both  the  particle  code  and 
the  prescribed  current  code.  Prescribed  ground  currents  without  a delay 
are  valid  in  the  ground  because  of  its  high  density  and  the  correspondingly 

3 

smaller  electron  lifetime  (shorter  by  a factor  of  10  ).  Conductivity 
enhancement  was  not  used  in  these  calculations. 


The  effect  of  ground  drivers  is  shown  in  the  comparison  of  Figure 
30.  The  calculations  shown  here  are  for  a fat  pole  with  enhancement  (particle 
sources) . Unfortunately,  there  was  an  error  in  the  enhancement  such  that  it 
was  a factor  of  1.5  times  greater  in  the  calculation  without  drivers  than 
in  the  one  with  drivers.  Our  previous  results  show  that  this  fact  will  not 
substantially  alter  the  comparison  however.  The  currents  seen  at  +30  cm, 
which  were  not  altered  by  the  presence  of  enhancement  at  all,  show  a 
slight  decrease  when  drivers  are  removed.  The  current  calculated  at  +70 
cm  showed  no  change  at  all,  and  the  comparison  is  not  shown.  In  the  ground, 
we  see  that  the  large  positive  peak  was  due  entirely  to  the  ground  drivers. 
Without  them,  we  have  a dominant  negative  peak  followed  by  a gradual  flow  of 
electrons  back  up  to  the  surface. 
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Comparison  between  pole  currents  predicted  with  and  without 
ground  currents.  Enhancement  in  problem  without  currents  Is 
1.5  times  greater  than  in  problem  with  ground  currents,  ex- 
plaining difference  in  late  time  pole  currents  below  surface. 
There  was  no  difference  In  pole  current  observed  at  +70  cm. 
This  calculation  used  particle  sources  in  the  air. 


The  part  of  the  current  identified  with  the  ground  drivers  is  very 
insensitive  to  the  pole  size,  in  the  two  cases  considered  here.  This  is 
true  whether  or  not  ground  conductivity  enhancement  is  considered.  A 
comparison  of  Figures  26  and  27  shows  that  this  is  true  without  enhancement 
and  the  comparison  of  Figure  31  shows  that  this  is  true  with  enhancement. 

This  insensitivity  is  consistent  with  the  idea  that  current  is  being  drawn 
from  an  area  of  ground  with  a radius  on  the  order  of  a skin  depth.  The 
differences  in  current  would  be  due  to  the  relative  areas  occupied  by  the  pole 
and  the  radial  electric  fields  developed  around  the  pole.  The  current  seen 
at  -45  cm  is  not  due  to  any  Compton  charge  generated  at  the  same  depth, 
but  is  the  reaction  to  fields  generated  by  the  radial  ground  conduction 
currents . 

As  mentioned  above,  a pole  transparent  to  gamma  rays  was  used  in 
these  calculations.  Similar  calculations  were  made  for  an  opaque  pole, 
but  no  significant  differences  were  seen.  This  may  be  due  in  part  to  the 
pole  size.  However,  there  is  good  reason  to  believe  that  the  actual  pole 
response  is  insensitive  to  pole  opacity  at  these  incidence  angles  and 
gamma  fluxes.  However,  the  fields  in  space  aiound  the  pole  are  sensitive 
as  was  seen  in  the  Ex  study  (Section  4.2). 

The  pole  response  due  to  the  direct  interaction  of  the  gamma  rays 
with  the  pole  obviously  become  reduced  as  the  incidence  angle  increases  and 
the  horizontal  component  of  the  flux  decreases.  There  is  an  additional 
reason  for  direct  interaction  effects  to  be  reduced  as  the  incidence  angle 
increases  in  the  presence  of  a large  gamma  pulse  with  its  correspondingly 
high  air  conductivity.  A gamma  wave  front  striking  the  pole  produces  either 
a new  increase  or  decrease  of  charge  on  it  depending  on  several  factors  which 
include  the  opacity  of  the  pole.  The  charge  can  be  neutralized  by  currents 
running  on  the  pole  and  by  conduction  currents  in  the  £ If  the  neutraliza- 
tion is  produced  by  currents  running  along  the  length  of  the  pole,  the  signature 
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Figure  31.  Comparison  of  currents  calculated  on  the  fat  and 

narrow  poles  by  DAVID  (particle  sources).  Both  cal 
culations  include  ground  conductivity  enhancement. 
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of  the  net  current  will  indicate  the  pole  opacity.  If  air  conduction  cur- 
rents dominate,  they  can  reduce  any  longitudinal  pole  currents  due  to  the 
direct  interaction  mechanism.  The  fluxes  used  in  this  study  produce  a high 
conductivity  very  early  in  the  pulse.  The  direct  interaction  signal  is  wiped 
out  before  it  has  a chance  to  be  seen  (actually,  there  was  a factor  of  two 
increase  in  the  tiny  scattered  field  response  of  the  fat  pole  at  +30  cm) . 

When  the  angle  of  incidence  is  small,  the  direct  interaction  occurs 
simultaneously  along  the  pole.  There  is  no  resulting  longitudinal  current 
•at  a given  point  until  such  time  as  a pulse  could  travel  from  some  discontinuity 
such  as  the  end  of  the  pole  or  the  air/ground  interface.  If  the  conductivity 
builds  up  in  this  time  frame,  the  propagated  longitudinal  current  pulse  will 
never  reach  the  observer.  This  shallow  incidence  angle  effect  was  probably 
another  factor  in  the  loss  of  an  opacity  signature  in  these  calculations. 

For  a given  peak  flux,  there  should  be  a incidence  angle  which  maximizes 
direct  interaction  effects,  as  seen  in  the  longitudinal  current. 


Figure  32  is  our  estimate  of  the  net  longitudinal  currents  that 
would  be  seen  at  +30  cm  and  -45  cm  on  a 10  cm  diameter  pole.  It  is  obtained 
from  the  particle  calculation  for  the  narrow  pole,  considering  both  ground 
drivers  and  ground  conductivity  enhancement.  The  following  types  of  scaling 
were  performed:  (1)  the  +30  cm  current  was  scaled  by  the  ratio  of  the  pole 

circumferences  (0.523),  (2)  the  negative  scattered  field  response  at  -45  cm 
was  scaled  by  the  same  factor,  and  (3)  the  positive  (ground  driver)  response 
at  -45  cm  was  not  scaled  at  all. 


Experimental  results  indicate  a peak  response  in  the  air  (+20  cm) 
of  2000  amps  and  in  the  ground  (-50  cm)  of  1900  amps,  with  the  peak  of  the 
ground  signal  occurring  about  10  ns  earlier  than  the  air  peak.  Our  pre- 
dictions (at  somewhat  different  positions)  indicate  peaks  of  1800  amp  and 
2500  amp  respectively  with  the  difference  in  peak  times  being  about  15  ns. 

Given  the  physical  and  numerical  uncertainties  involved,  this  agreement  must 
be  considered  fortuitous,  especially  in  the  ground.  The  currents  in  the  air, 
at  this  incidence  angle,  are  relatively  insensitive  and  should  have  been  easier 
to  estimate. 
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Figures  33  through  35  show  the  unsealed  currents  running  on  each 
side  of  the  narrow  pole.  These  can  be  quite  different  and  it  is  good  to 
remember  this  when  we  speak  of  the  net  longitudinal  current.  There  is  also 
a strong  current  running  around  the  side  of  the  pole  in  the  z-direction. 

We  did  not  calculate  it  explicitly,  but  its  existence  can  be  seen  from  the 
fact  that  H is  the  same  magnitude  as  H at  a position  Ax/2  away  from  side  1. 

y ^ 

To  obtain  the  net  longitudinal  current,  add  the  components  labeled  "front,” 
"back,"  and  "side"  with  the  indicated  sign  and  multiply  by  two.  The  values 
plotted  are  for  the  half  of  the  pole  on  one  side  of  the  symmetry  plane. 

When  comparing  the  current  components,  remember  that  the  width  of  pole 
over  which  the  side  current  density  is  integrated  is  four  times  greater 
than  the  width  of  pole  that  the  front  and  back  currents  were  integrated. 


36d  compare  the  spatial  variations  of  E and  H in 

y * 


Figures  36a 

the  z-direction  with  and  without  the  small  pole.  The  comparisons  are  at 
10.2  and  20.2  nsec.  The  variation  is  along  a line  30  cm  above  the  ground 
and  in  the  symmetry  plane.  Particle  sources  are  used.  The  pole  calculation 
is  the  same  computer  run  that  generated  the  currents  shown  in  Figures  32-35. 

Our  final  parameter  study  is  intended  to  show  the  importance  of 
the  scattered  field  response  for  small  incidence  angles  and  low  fluxes.  We 
do  not  claim  that  the  predictions  of  the  3-D  code  are  correct  here,  since 
it  does  not  calculate  the  propagated  signal  well,  but  it  should  certainly 
show  the  proper  trends.  The  analytic  source  code  (DAVEJR)  was  used  for 
these  studies. 
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Figure  37  shows  the  narrow  pole  currents  calculated  with  a neak 
13 

dose  of  10  rad/sec  and  an  incidence  angle  of  20°.  Ground  drivers  and 
conductivity  enhancement  are  included.  In  Figure  38,  we  show  the  early  time 
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Figure  34 


. Components  of  current  running  along  sides  of  the  "narrow" 
pole  at  +70  cm. 
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Figure  36a.  Variation  of  magnetic  field  in  z-direction  with  pole  (WP)  and 
without  pole  (WOP)  along  a line  at  +30  cm  height  on  symmetry 
plane.  Particle  sources,  T * 10.2  ns.  1 gauss  - 30,000  v/m. 


Figure  36b.  Variation  of  vertical  electric  field  in  z-direction  with  pole 
(WP)  and  without  pole  (WOP)  along  a line  at  +30  cm  height  on 
symmetry  plane.  Particle  sources,  T * 10.2  ns. 
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Figure  36c.  Variation  of  magnetic  field  in  z-direction  with  pole  (WP)  and 
without  pole  (WOP)  along  a line  at  +30  cm  height  on  symmetry 
plane.  Particle  sources,  T » 20.2  ns. 
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Figure  36d.  Variation  of  vertical  electric  field  in  z-direction  with  pole 
(WP)  and  without  pole  (WOP)  along  a line  at  +30  cm  height  on 
symmetry  plane.  Particle  sources,  T = 20.2  ns. 
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portion  of  the  pulse  when  the  incidence  angle  is  reduced  to  36.  Note  the 
prominance  of  the  negative  scattered  field  response,  even  at  these  dose 
rates.  In  Figures  39  and  40  we  maintain  the  same  angle  (38)  but  reduce 
the  flux  to  10  rad/sec  and  4 • 10  rad/sec  respectively.  Even  though 
the  calculations  were  only  taken  out  to  30  nsec,  one  can  see  that  the 
scattered  field  response  has  become  totally  dominant  (except  at  +70  cm). 
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SECTION  5 

CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  described  a three-dimensional  close-in  coupling  code  con- 
structed for  the  purpose  of  estimating  the  voltages  and  currents  induced 
on  arbitrarily  shaped  objects.  The  code  comes  in  two  versions:  a particle 

following  self-consistent  version  (DAVID),  and  a simpler  one  which  uses 
prescribed  currents  in  the  air  (DAVEJR).  The  two  versions  are  identical 
except  for  the  source  routines.  Since  DAVEJR  does  not  store  particle  infor- 
mation or  solve  the  equations  of  motion,  it  is  much  faster  and  can  be 
optimized  to  handle  a much  larger  grid.  It  should  be  an  extremely  useful 
tool  in  the  study  of  the  tactical  problem. 

The  philosophy  used  in  writing  DAVID  and  DAVEJR  emphasized 
simplicity  and  usability.  The  codes  could  have  been  made  more  sophisticated, 
both  in  terms  of  the  physics  and  the  numerics.  Instead,  an  attempt  was 
made  to  produce  a reliable  code  with  well  defined  and  tested  limits  of 
applicability;  one  which  could  easily  be  upgraded  or  modified  by  future 
users.  This  goal  seems  to  have  been  accomplished.  The  codes  give  the 
right  answers,  within  the  limits  built  into  them.  The  coding  is  relatively 
easy  to  follow  and  other  users  should  have  no  great  difficulty  in  adapting 
DAVID  to  their  needs. 

DAVID  was  tested  against  the  GLANC  1-D  close-in  environment  code. 
The  results  are  reported  in  Section  4.2.  GLANC  was  configured  to  match 
the  DAVID  calculation  in  terms  of  grid  and  time  step  size.  In  general,  the 
agreement  was  quite  good.  The  disagreements  were  explained.  One  result 


of  this  comparison  was  to  increase  the  working  volume  before  we  attempted  the 
pole  current  calculations  shown  in  Section  4.3.  The  results  of  the  pole 
current  parameter  study  were  explained  theoretically  and  gave  insight  into 
how  one  might  construct  models  or  use  simpler  techniques  to  calculate  such 
problems  in  the  future. 

Recommendations  for  future  work  fall  into  two  categories:  (1)  types 

of  problems  to  be  attacked  and  (2)  code  improvements.  The  types  of  improve- 
ments will  depend,  in  part,  on  what  problems  are  solved. 

In  many  ways,  3-D  close-in  coupling  codes  are  very  powerful. 

They  are  also  very  limited.  They  include  a great  deal  of  physics,  but  lack 
spatial  resolution,  are  limited  to  fairly  high  conducting  regions,  and  are 
expensive  to  use.  The  3-D  code  is  but  one  weapon  in  the  arsenal  we  have  at 
our  disposal  to  attack  the  close-in  coupling  problem.  The  other  weapons 
include  analytic  calculations,  lumped  and  distributed  parameter  models, 
and  2-D  codes.  DAVID  can  best  be  utilized  in  (1)  calculating  currents  and 
voltages  on  arbitrarily  shaped  objects,  within  the  constraints  of  the  code, 

(2)  isolating  the  important  parts  of  the  physics  so  that  other  (simpler) 
techniques  can  be  developed,  and  (3)  aiding  in  the  development  of  methods 
for  building  lumped  and  distributed  parameter  models  of  complicated  systems 
and  their  interaction  with  the  environment.  All  too  often,  codes  of  this 
nature  are  considered  ends  in  themselves.  They  become  crutches  and  an 
excuse  to  stop  thinking.  The  state-of-the-art  is  not  nearly  advanced  enough 
to  build  black  box  type  codes. 

In  the  immediate  future,  DAVID  (and  DAVEJR)  should  be  used  to  continue 
the  theoretical  study  of  the  pole  problem  begun  in  Section  4-3.  This  should  lead 
into  the  low  flux  "tactical"  problem.  It  would  be  very  useful  to  introduce 
variously  shaped  objects  and  to  look  at  the  effect  one  nearby  object  has  on 
another.  One  could  look  at  two  poles  in  different  relative  positions  and 
even  connect  them  to  form  loops  with  a magnetic  field  response. 


The  pole  calculation  pointed  out  an  important  factor  that  must 
be  considered  when  going  to  tactical  problems.  That  factor  is  the 
importance  of  the  scattered  field  response  of  the  pole  to  the  propagated 
vertical  electric  fiel 1 and  the  vertical  component  of  the  local  radial  (from 
the  burst)  electric  field.  This  response  is  very  important  below  the  ground 
even  with  high  fluxes,  and  in  the  air  when  the  angle  of  gamma  incidence  is 
small  and/or  the  flux  is  small.  The  three-dimensional  codes  do  not  calculate 
this  part  of  the  response  well  because  a long  distance  is  required  to  build 
up  the  proper  propagated  field*.  One  could  put  in  a more  complicated 
boundary  condition  using  a 2-D  environment  code  calculation,  but  a wiser 
move  would  probably  be  to  do  the  response  problem  in  two  steps.  The  early 
time,  low  conductivity,  part  of  the  response  can  be  done  using  other  methods, 
while  the  later,  high  conductivity  part  of  the  calculation  can  be  done  using 
a full  3-D  calculation. 

The  low  conductivity  scattered  field  response  could  be  obtained  in 
one  of  several  ways.  Analytic  calculations  might  fc-  appropriate.  On  the 
other  hand,  if  the  pole  is  thin  compared  to  interesting  wavelengths,  a 2-D 
code  might  be  used  in  which  minus  the  incident  vertical  electric  field  is 
used  as  the  pole  boundary  condition  and  the  pole  is  allowed  to  radiate  into 
a time  dependent  conductivity,  which  is  obtained  from  an  environment  code 
calculation.  Codes  of  this  nature  already  exist  and  are  easily  constructed 
in  any  case.  For  more  complicated  structures,  a 3-D  code  can  be  easily 
and  cheaply  constructed  using  the  same  principle**.  Such  a 3-D  code  could 


* Cylindrical  geometry  codes  do  even  more  poorly  because  the  grid  points 
do  not  line  up  along  the  wave  front  and  because  the  cells  become  so 
large  near  the  outer  boundary. 

**  The  appropriate  parts  of  DAVID  or  DAVEJR  could  be  used  to  greatly  reduce 
the  effort  required  to  build  a new  code. 


have  a very  large  grid  and  would  be  useful  for  doing  scattering  problems 
even  in  the  absence  of  any  conductivity.  DAVIDS  field  equations  are  dif- 
ferenced in  such  a way  as  to  allow  zero  conductivity  even  now.  They  could 
be  written  to  remove  the  conductivity  term  and  thus  be  even  faster. 

General  improvements  which  could  be  made  in  DAVID  and  DAVEJR  at 
this  time  involve  the  outer  boundary  condition  and  the  grid  spacing.  It 
would  be  useful  to  provide  layers  of  smaller  sized  cells  above  and  below 
the  ground  plane  and  in  the  immediate  vicinity  of  the  object.  Large 
cells  can  be  used  away  from  those  two  regions.  The  outer  boundary  condition 
involves  both  the  particles  and  the  electromagnetic  field.  A 1-D  environ- 
ment calculation  could  be  used  to  supply  the  fields  and  a particle  distribu- 
tion. Only  the  field  boundary  condition  would  be  required  in  DAVEJR,  where 
prescribed  sources  are  used.  The  present  boundary  condition,  which  uses 
the  radial  electric  field,  is  described  in  Section  3.2.  As  mentioned  above, 
one  could  use  a 2-D  environment  calculation  at  tne  front  boundary  to  include 
the  propagated  field,  but  then  one  would  also  have  to  worry  about  the 
scattered  field  from  the  object.  Other  methods  should  be  pursued  first. 

A carefully  considered  treatment  of  electrons  emitted  from  the  pole 
and  backscattered  from  the  soil  should  be  implemented.  It  is  very  important 
that  such  a treatment  be  consistent  with  the  level  of  the  physics  already 
included  so  that  the  proper  ratio  of  currents  is  maintained. 

The  DAVII)  concept  should  also  be  used  to  calculate  system  responses 
in  the  high-altitude  burst  close-in  region,  or,  to  some  degree,  in  the 
investigation  of  3-D  effects  in  the  SGEMP  problem.  There  are  several  problems 
inherent  in  the  last  application  which  involve  the  definition  of  currents 
near  the  satellite  and  the  relatively  small  struts  which  support  the  solar 
panels.  It  would  be  very  useful  to  build  a simple  code  which  used  non- 
self-consistent  prescribed  currents  or  prescribed  currents  with  an  approxima- 
tion for  self-cor.sistency. 
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A listing  of  DAVID,  DAVEJR,  and  the  time  waveform  output  code 
DAVEOUT  can  be  obtained  by  qualified  users  through  the  Defense  Nuclear 
Agency. 
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APPENDIX  A 

AN  APPROXIMATION  FOR  INCLUDING  MAGNETIC  TURNING 
EFFECTS  IN  PRESCRIBED  ELECTRON  CURRENTS 


Efforts  are  now  being  made  to  perform  three-dimensional  close-in 
EMP  calculations.  In  the  regions  of  interest,  the  electric  and  magnetic 
fields  are  great  enough  to  influence  the  motion  of  the  Compton  electrons, 
which  originally  created  them.  Self-consistent  particle  following  schemes 
are  the  best  way  to  calculate  these  driving  currents.  However,  particle 
calculations  require  considerable  computer  storage  and  time.  This  time 
and  storage  can  often  be  put  to  better  use,  e.g.,  proving  a more  accurate 
system  model  or  generating  more  detailed  parameter  studies.  A prescribed 
current  calculation  can  be  used  to  generate  a current  description  at  low 
altitudes.  Using  this  technique,  the  current  is  given  by  the  gamma  flux 
times  a constant  (Reference  A-l).  The  same  is  true  of  the  ionization  rate. 
It  seems  reasonable  to  expect  that  an  approximation  should  be  available 
which  would  give  a first  order  correction  to  the  prescribed  current  to 
account  for  the  effects  of  magnetic  turning  and  the  electric  field  drag 
force.  One  such  technique  (Reference  A-2),  which  we  will  call  the  "DX-DY" 
approximation,  already  exists.  It  is  a good  approximation.  There  is  a 
simpler  way  to  obtain  useful  estimates,  however,  a method  which  can  easily 
be  made  more  accurate  after  a few  parameter  studies  using  particle  pushing 
codes  have  been  made. 


The  most  important  self-consistent  effect  for  fields  near  a sur- 
face at  sea  level  is  the  magnetic  turning  effect.  Compton  electrons 
traveling  initially  parallel  to  the  ground  are  turned  upward,  driving  a 
significant  electric  field  which  couples  rather  efficiently  into  objects  poking 
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up  out  of  the  earth.  The  normal  drag  forces  more-or-less  mask  the  electric 
field  drag.  Therefore,  electric  field  effects  will  be  ignored  here,  al- 
though they  could  be  included  with  no  difficulty.  Actually,  only  one  mag- 
netic field  component  will  be  considered:  the  one  parallel  to  the  ground 

and  the  incident  wave  front  (the  x-direction  in  Figure  A-l).  This  is  the 
component  which  would  exist  in  the  absence  of  an  object  and  which  is  normal 
to  the  path  of  the  incoming  electrons.  The  current  running  on  a vertical 
pole  will  generate  an  azimuthal  field  which  adds  to  the  ambient  field  on 
one  side  and  subtracts  from  it  on  the  other. 

We  start  with  the  equation  of  motion  of  an  electron  in  a magnetic 
field,  w’th  collisions: 


•dT+vPz  = Vy 


dt  + vPy  = • "'"n\  • 


(A-2) 


where 


mVl  + (P/mc)2 


(A-3) 


V*  • 


. . x io8o  (v°-» 

312  E (E  +0.6)  ’ 

» e e 

2 2 

P = P + P , magnitude  of  momentum 
y z 

E = electron  kinetic  energy  (MeV) 

-19 

e = electron  charge  (1.602  x 10  coulomb) 

3 

p = air  density  (1.23  kg/m  at  sea  level) 


(A-5) 


m = electron  mass  (9.108  * lO-15  kg) 
Bx  = magnetic  field  intensity  (Tesla) 
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Figure  A-1.  Problem  geometry.  Gamma  wave  front  Is  parallel  to 
x-axis.  Wave  normal  makes  angle  e,  measured  from 
the  horizontal. 

The  quantities  v and  will  be  treated  as  constants.  Representative  values 
of  P and  E must  be  used  to  calculate  them  directly.  They  will  actually  be 

v 

found  from  comparison  with  more  accurate  calculations  as  a function  of  gamma 
energy. 


We  define  the  complex  quantity 

P “ P2  + iPy  • (A-6) 

Equations  A-1  and  A-2  combine  to  yield 
dP 

^ + (V  - iu>H)P  = 0 , (A-7) 


which  has  the  solution 


where  y = v - iu^  and  PQ  = PQz  + iPQy,  the  initial  momenta.  The  average 
value  of  P over  the  time  interval  At  is 


PA  * y&  V - e'YAt)  ' CA-9) 

We  will  be  interested  in  how  the  avera-e  momentum  with  turning  compares  to 
the  average  momentum  without  turning,  since  this  will  indicate  how  the  cur- 
rents compare.  First,  assume  that  vAt  is  large,  which  rids  us  of  the  damped 
oscillation  contained  within  the  parentheses.  Then,  the  real  and  imaginary 
parts  of  the  modified  equation  yield 


U>HPyO 


vPzO  - V 
(v2+u£)At 


(A- 10) 


(v2+u^)At 


CA— 11) 


where  P&z  and  are  the  average  momentum  components.  Without  turning, 
these  components  would  be 


P. 

az  vAt  ’ 


p,  . !°jr 

ay  vAt 


Then,  in  terms  of  the  non-turning  momentum, 


(A- 12) 


(A-13) 
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(A-15) 
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We  now  assume  for  curve-fitting  purposes  that  the  primed  and  un- 
primed momenta  correspond  to  the  standard  prescribed  current  and  tne  modified 
current.  Let  J be  the  prescribed  current  that  would  be  calculated  without  a 
magnetic  field,  and  let  Jsc  be  the  current  calculated  with  a magnetic  field. 
Then,  we  propose  that  they  are  related  in  the  following  manner: 


- * 


(A- 16) 


J +■ 


1 -4-..- 


(A-17) 


Note  that  the  self-consistent  current  depends  upon  only  the  ratio 
(u^/v),  which  can  be  written  as  (BX/BA),  i.e., 

T E 17  • tA-18^ 

A 

where  is  a function  of  the  electron  energy  and  is  proportional  to  the 
air  density.  It  has  units  of  magnetic  field  intensity.  It  will  be  fitted 
as  a function  of  the  gamma  ray  energy. 

To  simplify  the  equations,  rotate  the  coordinate  axis  such  that 
the  z-axis  is  parallel  to  the  direction  of  the  gamma  rays.  We  can  then 
refer  to  the  original  prescribed  current  as  J (it  being  parallel  to  the 
z-axis).  The  component  of  the  self-consistent  current  in  this  direction 
will  be  called  JD  and  the  component  turned  normal  to  J will  be  called  JM. 
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(A- 20) 


(A-21) 


These  ratios  were  calculated  as  a function  of  gamma  ray  energy 
(assuming  that  the  Compton  electron's  initial  energy  was  half  the  gamma 
energy)  and  Bx<  The  initial  momentum  was  calculated  through 


PQ  * (me) 


[(■•5)’ 


1/2 


(A-22) 


The  drag  and  turning  constants  (lumped  into  BA)  were  calculated  using  the 
initial  value  of  momentum  and  energy.  This  leads  to  smaller  values  of  v 
(and  B^)  than  are  realistic,  but  the  ratios  will  still  be  reasonable.  The 
ratios  are  shown  in  Figures  A-2,  A-3,  and  A-4  as  a function  of  Bx  and  the 
gamma  ray  energy.  They  are  compared  with  the  same  ratios  computed  from  the 
DX-DY  data  published  in  Reference  A-2*.  The  agreement  is  within  a factor 
of  two  over  most  of  the  parameter  range.  The  new  results  are  much  smoother, 
however.  This  is  not  an  indication  that  they  reflect  physical  reality  any 
better.  Both  methods  become  increasingly  inaccurate  as  the  magnetic  field 
increases.  The  smoother  curves  generated  by  the  hew  technique  aid  in  the 
numerical  integration  of  Maxwell's  equations,  however. 


In  order  to  improve  the  approximation,  we  decided  to  calculate 
using  Longley’s  (Reference  A-2)  DX-DY  data,  rather  than  trying  to  choose 
an  average  value  of  electron  momentum  or  some  other  such  thing.  The  DX-DY 
predictions  are  most  accurate  at  low  magnetic  fields.  Since  only  one 


* Reference  A-2  contains  data  for  a wide  range  of  gamma  energies,  electric 
fields,  and  magnetic  fields. 
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New  Approximation 
Using  Initial 
Electric  Momentum 
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Figure  A-2.  Ratio  Op/J  as  a function  of  magnetic  field  intensity  for 
three  gamma  ray  energies  (1.0,  1.5,  2.0  MeV).  Curves  are 
shown  for  the  DX-DY  approximation,  the  new  approximation, 
and  the  new  approximation  fitted  to  the  DX-DY  predictions 
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Figure  A-4.  Ratio  0m/0p  as  a function  of  magnetic  field  intensity  for 
three  gamma  energies  (1.0,  1.5,  2.0  MeV).  Curves  are  shown 
for  the  DX-DY  approximation,  the  new  approximation,  and  the 
new  approximation  fitted  to  the  DX-DY  predictions  (1  MeV). 


parameter,  B^,  is  used  in  the  new  approximation,  it  was  easy  to  compute  it 
from  the  ratios  published  in  Reference  A-2  for  zero  electric  field.  BA 
is  given  by  Equation  A-21.  The  ratio  corresponding  to  J^/Jp  in  Reference  A- 
2 is  DY/DX.  Table  A-l  shows  the  fitted  values  of  (in  gauss)  as  a func- 
tion of  incident  gamma  energy,  E^.  E^  is  between  0.5  MeV  and  6.0  MeV. 
Fortunately,  a rather  simple  curve-fit  of  this  data  exists.  It  is 

g 

BA(V  " in(E°/E0')'  * 0.5  < Ey  < 6 MeV  , (A-23) 

where 

Bq  = 44.27  Pr  (gauss) 

EQ  = 0.3322  (MeV) 

Pr  = air  density,  relative  to  sea  level 
E^  = gamma  energy  (MeV) 

An  example  of  the  new  fitted  calculation  is  also  shown  in  Figures 
A-2  through  A-4.  The  example  is  that  of  a 1 MeV  gamma  ray. 


In  DAVEJR,  and  Jp  are  calculated  using  Equations  A-19  and 
A-20  and  a coordinate  system  rotation,  when  the  parameter  ISELF  is  set  equal 
to  1.  When  ISELF  « 0,  no  self-consistent  effects  are  considered.  There  is 
also  a smearing  of  the  H-field  (HX),  which  is  explained  later. 


The  prescribed  current  formulation  assumes  that  the  electron  cur- 
rent is  in  equilibrium  with  the  gamma  flux,  i.e.,  the  average  Compton 
collision  frequency  is  high  enough  to  erase  any  of  the  electrons  memory  of 
previous  changes  in  the  gamma  flux  time  history  over  its  lifetime.  With 
turning  being  added  to  the  current  description  in  the  way  that  it  was  done 
here,  that  assumption  is  modified  by  the  assumption  that  the  lifetime  of  the 
electron  is  short  compared  to  the  time  that  the  magnetic  field  changes 
significantly. 
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Table  A-l.  Ba  as  a function  of  gamma  energy  (Ey). 

These  values  are  calculated  using  tne 
OX-DY  data  published  In  Reference  A-2. 


Ey(MeV) 

BA(gauss) 

0.5 

84.00 

1.0 

40.17 

1.5 

29.47 

2.0 

24.66 

3.0 

20.15 

4.0 

17.97 

5.0 

16.65 

6.0 

15.80 

A word  of  caution  must  be  expressed.  When  calculating  currents 
near  a surface,  one  must  remember  that  the  surface  will  prevent  the  complete 
turning  of  the  electrons  and  hence  will  reduce  the  influence  of  the  magnetic 
field  in  changing  the  prescribed  current.  The  surface  effect  will  increase 
for  increasing  incidence  angles,  but  will  decrease  for  increasing  magnetic 
field  strengths,  since  the  turning  radius  of  the  electron  decreases  and  one 
can  get  closer  to  the  surface  before  electrons  collide  (see  Figure  A-5). 
Obviously,  the  range  of  the  electron  is  also  important  in  determining  the 
height  at  which  the  ground  influence  will  be  felt.  To  a lesser  degree,  the 
angular  distribution  of  the  Compton  electrons  is  of  interest.  At  near 
grasing  angles,  the  fact  that  Compton  electrons  are  ejected  in  a cone  about 
the  direction  of  gamma  incidence  contributes  to  a small  net  vertical 
current.  None  of  these  effects  have  been  considered  here.  A relatively 
simple  height  dependence  can  probably  be  built  into  the  "constant"  BA  to 
obtain  an  estimate  of  the  importance  of  the  surface  interaction  (which  also 
includes  backscattered  electrons,  another  effect  that  we  have  ignored). 


i 
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Figure  A-5.  Effect  of  a surface  on  the  current 
components  generated  by  magnetic 
turning. 


A second  problem  occurs  with  the  use  of  modified  prescribed 
sources.  The  approximations  assume  that  the  electron  sees  a constant 
magnetic  field  (and  electric  field,  when  used)  over  its  lifetime.  In  sea 
level  air,  that  time  is  about  1 shake.  Thus,  an  electron  bom  near  the 
beginning  of  the  pulse  can  experience  a wide  range  of  field  strengths  over 
its  lifetime. 

In  DAVEJR,  we  multiply  the  HX  field  by  a time-dependent  quantity 
which  reduces  the  effective  H-field  during  the  rise  of  the  pulse.  At  late 
times,  when  the  magnetic  field  is  presumably  varying  slowly,  the  function 
approaches  unity  and  the  prescribed  current  then  reacts  to  the  instantaneous 
field.  The  function  is  called  HSMR  (H-smear)  and  is 


HSMR  = 


1 - EXP(-ATSMR) 
ATSMR 


CA-24) 


where 


ATSMR  = y [1  - TANH (ALPHA  • TSMR)]  (ALPHA  • DTSMR)  + 10-4 
ALPHA  = rise  e-fold  rate  of  pulse 

DTSMR  = time  over  which  the  magnetic  field  is  to  be  smeared 
TSMR  = TS-TPK 


TPK  = time  at  which  gamma  pulse  peaks 

TS  = local  time  at  which  source  is  being  calculated. 

This  is  not  a true  "smearing"  of  the  magnetic  field,  which  would  involve 
a convolution  with  a memory  function  or  a weighted  average.  Such  techniques 
involve  storing  past  values  or  functions  of  the  field.  With  a three- 
dimensional  code,  this  requires  significant  storage.  In  order  to  avoid  the 
use  of  this  storage,  the  method  using  Equation  A-24  was  devised. 
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